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1  INTRODUCTION 


The  design  of  practical  combustors  has  been  and  still  is  for  the  most  part  an  art[l]. 
Traditional  design  procedures  are  based  on  “cut-and-try”  methodologies  and  are 
guided  by  experienced  combustor  designers.  The  complexity  of  the  fluid  dynam¬ 
ics  and  chemistry  has  defied  rigorous  detailed  numeiical  simulation.  In  practical 
applications  even  axisyrninctric  dump  combustors  as  shown  in  Figure  1  exhibit  very 
complex  flow  fields.  They  usually  have  three-dimensional  characteristics  due  to  swirler 
vanes  or  injection  of  fuel  at  discrete  locations.  The  flow  may  have  several  gaseous 
components,  liquid  droplets,  and  solid  particles  that  physically  and  chemically  react. 
Thermal  radiation  may  be  important  with  chemical  reaction.  Swirl  flows  can  create 
a  central  toroidal  recirculation  zone  (CTRZ)  in  addition  to  the  corner  recirculation 
zone  (CRZ)  and  hub  recirculation  zone  (HRZ).  In  this  study  the  effects  of  a  closely 
coupled  CD  nozzle  also  complicates  the  flow  field.  Experimental  studies[2,  3]  show 
that  the  existence,  shape,  and  strength  of  the  CTRZ  is  strongly  influenced  by  details 
of  the  combustor  geometry,  upstream  flow  conditions,  downstream  flow  conditions, 
combustion,  and  amount  of  swirl. 

The  sudden  expansion  dump  combustor  illustrated  in  Figure  1  is  favored  for  use  in 
ramjet  propulsion  systems.  The  dump  combustor  is  simple  and  can  be  designed  to  be 
short  and  compact,  especially  if  swirl  is  used  to  stabilize  the  combustion  flame.  Swirl 
is  induced  into  the  flow  entering  the  combustor  by  use  of  turning  vanes  or  tangential 
injection  of  air.  The  amount  of  swirl  is  measured  by  the  swirl  number  (S)  which  is 
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Figure  1:  Schematic  Diagram  of  a  Dump  Combustor  Flow  Field 
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The  swirl  number  is  approximately  equal  to  two-thirds  of  the  tangent  of  the  swirl 
vane  angle  (^swirj)-  Swirl  has  a  strong  effect  upon  the  combustor  flow  ficld[3].  At 
low  degrees  of  swirl  ( S  <  0.4)  the  width  of  the  air  jet  is  increased  over  that  without 
swirl.  At  larger  swirl  numbers  the  mixing  rate  of  the  jet  increases.  At  swirl  num¬ 
bers  of  approximately  0.6  and  greater  a  CTRZ  forms.  Formation  of  the  CTRZ  is 
usually  described  as  “vortex  breakdown”  or  a  transition  from  supercritical  swirl  flow 
to  subcritical  swirl  flow[3,  4j.  The  CTRZ  can  be  designed  to  be  very  stable  and  to 
provide  flame  holding  capability.  The  reverse  flow  mass  flowrate  in  the  CTRZ  can 
be  larger  than  the  mass  flowrate  in  the  primary  air  jet.  The  enhanced  mixing  rate  in 
the  CTRZ  improves  combustion  and  thus  allows  the  length  of  the  dump  combustor 
to  be  shortened.  This  saves  weight  and  space  which  can  be  important  in  designing 
propulsion  systems.  In  addition,  using  the  CTRZ  for  flameholding  has  advantages 
over  bluff  body  flame  holders  and  simple  dump  combustors  (without  swirl)  in  that 
the  hot  region  of  the  flow  is  not  near  a  wall  thereby  reducing  cooling  requirements. 

Making  dump  combustors  shorter  increases  the  influence  of  the  exhaust  nozzle  upon 
the  combustor  flow.  Particularly  for  high  swirl  numbers  the  exhaust  nozzle  can  have 
a  substantial  effect  on  the  CTRZ.  See  for  example  the  experimental  work  of  Escudier 
and  Keller[4]  showing  substantial  changes  in  axial  velocity  profiles  and  shape  of  the 
CTRZ  with  only  a  15%  reduction  in  exit  diameter  for  high  swirl  number  in  a  dump 
combustor.  Also  Lilley [5]  and  Yoon  and  Lilley[6]  report  effects  of  downstream  geom¬ 
etry  contractions  on  the  CTRZ.  Their  experimental  results  showed  that  the  length  of 
the  CTRZ  was  shortened  by  the  contraction  and  the  shape  of  the  CTRZ  substantially 
changed. 

Several  computer  codes  have  been  developed  for  simulating  the  fluid  flow  in  axisym- 
metric  dump  combustors.  A  review  of  some  of  these  is  given  by  Lilley [7] .  More  calcu¬ 
lations  have  been  presented  recently  by  Lilley [8],  Abujelala[9],  Ramos  and  Somer[10], 
and  Rhode  and  Stowers[ll]  and  others.  Nearly  all  of  these  recent  works  have  used 
a  variation  of  the  TEACH  codcli2]  or  a  code  using  a  similar  solution  strategy.  In 
these  codes  a  relaxation  procedure  is  used  to  solve  the  steady  time  averaged  Navier- 
Stokes  equations  with  models  for  turbulent  viscous  stresses.  The  strong  streamline 
curvature  and  recircul  „tion  regions  present  in  the  dump  combustor  requires  that  the 
full  set  of  elliptic  equations  be  solved  (as  opposed  to  a  spacial  marching  solution  of  a 
parabolized  subset  of  these  equations). 

typically  (as  doin’  iu  S  I’ARl’ICjS])  (lie  flowfield  variables  are  over  specified  at  the 
inflow  boundary  and  set  by  zero  axial  gradient  at  the  outflow  boundary.  Since  the 
equations  in  nearly  ail  of  these  above  mentioned  codes  are  essentially  incompressible 
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these  boundary  conditions  can  lead  to  correct  and  unique  physical  solutions.  Applica¬ 
tions  of  these  codes  to  compressible  flows,  however,  can  lead  to  nonphysical  solutions, 
since  the  boundary  conditions  are  not  consistent  with  method-of-characteristic  the¬ 
ory  and  the  absence  of  density  variations  (due  to  pressure  changes)  is  inconsistent 
with  transonic  flow.  Take  for  example  the  flow  through  a  choked  CD  (convergent- 
divergent)  nozzle.  The  specification  of  all  flow  field  variables  at  the  inflow  boundary 
determines  the  net  mass  flow  rate.  This  mass  flow  rate  will  not  necessarily  match  the 
critical  mass  flow  rate  through  the  CD  nozzle.  Information  in  the  form  of  acoustic 
waves  must  propagate  upstream  from  the  throat  of  the  CD  nozzle  and  adjust  the 
flow  conditions  at  the  inflow  boundary.  A  realistic  boundary  condition  would  be  to 
specify  total  pressure  and  total  temperature  at  the  inflow  boundary.  Acoustic  waves 
reaching  the  inflow  boundary  would  change  the  local  velocity  and  thus  the  local  static 
pressure  and  temperature  of  the  flow  until  the  physically  correct  steady-state  mass 
flow  rate  was  entering  the  combustor. 

Most  flow  analysis  codes  currently  being  used  for  computing  steady  turbulent  non¬ 
reacting  flows  in  axisymmetric  dump  combustors  are  based  upon  the  mathematical 
model  and  solution  procedure  of  the  TEACH  code.  The  specific  problem  addressed 
in  this  work  arises  from  the  fact  that  these  computer  codes  can  not  be  applied  to 
the  calculation  of  compressible  transonic  and  supersonic  flow  in  a  closely  coupled 
dump  combustor  and  exhaust  nozzle  configuration  due  to  inappropriate  boundary 
conditions  and  lack  of  compressibility  effects. 

Numerous  compressible  flow  Navier-Stokes  codes  do  exist  that  have  been  applied 
routinely  to  the  computation  of  flow  through  CD  nozzles  yielding  accurate  solutions. 
See  for  example  Peery  and  Forester[13]  and  Cline[14],  Codes  of  this  type  are  based 
upon  time-dependent,  methods  of  solving  the  unsteady  mass-averaged  Navier-Stokes 
equations.  Steady-state  solutions  are  obtained  by  integrating  the  time-dependent 
equations  until  the  solution  approaches  sufficiently  close  to  the  asymptotic  stead- 
state  limit.  Most  of  these  so  called  time-dependent  codes,  however,  do  not  include 
swirl  or  advanced  turbulence  models.  In  addition,  most  of  these  time-dependent  codes 
have  not  been  applied  routinely  to  low-speed  recirculating  flows.  The  reason  for  this 
is  that  most  of  these  codes  have  numerical  stability  constraints  that  limit  the  size  of 
the  time  step  to  a  CFL  number1  of  about  1  or  achieve  maximum  convergence  rate  (to 
steady  state)  at  a  CFL  number  of  about  1.  These  codes  have  used  various  explicit, 
semi-implicit,  and  ADI  implicit  solution  procedures.  To  reach  steady-state  in  a  low- 
speed  recirculating  flow  requires  thousands  if  not  tens  of  thousands  of  time  steps[15] 
using  these  time-dependent  codes. 

Recently  progress  has  been  made  in  improving  the  computational  efficiency  and  ac¬ 
curacy  of  solution  procedures  for  the  compressible  Navier-Stokes  equations.  Several 
investigators  including  MacConnackfl 6],  Thomas  and  Wa!ters[17],  Chakravarthy[18] 

'(’Ft  number  is  the  ratio  of  the  acoustic  length  scale  (u  4-  c, )6t  to  the  cell  size  Ax 


3 


have  used  line  Gauss-Seidel  (LGS)  relaxation  methods  with  second  or  third-order 
flux-split  differencing  on  the  convection  and  pressure  terms  and  second-order  central 
differencing  on  diffusion  terms.  The  methods  are  unconditionally  stable  on  model 
equations.  Convergence  rates  to  steady  state  have  beer  shown  to  be  much  less  sen¬ 
sitive  to  the  time  step  size  than  the  popular  approximately-factored  ( A !•')  implicit 
methods  ( e.g .  Beam  and  Warming(19]).  Once  past  starting  transients  MacCormack’s 
method  has  been  run  with  CFL  numbers  of  one  billion,  thus  reducing  to  a  modified 
Newton  iteration  scheme.  Convergence  to  steady  state  for  some  calculations[l6,  17] 
can  require  as  few  as  20  time  steps.  More  recently  Yoon [20]  has  developed  a  Lower- 
Upper  Symmetric  Gauss  Seidel  (LU-SGS)  solution  procedure  that  converges  rapidly 
and  is  computationally  efficient.  As  described  in  this  report  the  LU-SGS  scheme 
makes  it  possible  to  create  a  flow  analysis  that  has  less  numerical  dissipation  than 
the  LGS  methods  mentioned  above.  Both  the  LGS  and  the  LU-SGS  methods  were 
studied  and  compared  for  use  in  calculating  dump  combustor/CD  nozzle  flows. 

In  this  work  a  computer  program  was  developed  called  DUMPSTER.  It  solves  the 
Navier-Stokes  equations  for  axisymmetric  flows  in  dump  combustors  and  nozzles.  The 
flow  analysis  is  applicable  to  flows  with  low  subsonic,  transonic,  and  supersonic  re¬ 
gions.  It  optionally  solves  one  or  more  transport  equations  for  trace  chemical  species. 
Turbulence  is  modeled  by  either  the  Baldwing- Loinax[21]  algebraic  model  or  Launder 
and  Spalding’s  two-equation  k-e  model[22].  The  computational  grid  is  constructed 
from  stacking  one  or  more  computational  blocks  (or  zones)  to  simplify  the  specifica¬ 
tion  and  generation  of  the  computational  mesh  as  well  as  efficiently  utilize  all  of  the 
computational  cells.  DUMPSTER  was  applied  to  several  dump  combustor  flows  that 
also  had  experimental  measurements  to  which  the  calculated  flows  were  compared. 
Results  are  also  compared  to  calculations  performed  witli  a  different  computer  flow 
analysis  program. 
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2  MATHEMATICAL  MODEL 

2.1  Navier-Stokes  Equations 

The  DUMPSTER  code  solves  the  Reynolds-averaged  Navier-Stokes  equations,  trace 
species  transport  equations,  and  an  optional  two-equation  turbulence  model.  These 
equations  are  given  below  in  integral  form[23]  for  a  cylinderical  coordinate  system. 

mJfLudv+JJspndS+I!LNdv  =  0  (1) 

In  the  above  equation: 


t  is  the  time. 

V  is  the  volume  of  the  domain. 

U  is  the  solution  state  vector  u. 

S  is  the  surface  containing  the  volume. 

P  is  the  flux  vector. 

n  is  the  unit  normal  vector  pointing  outward  from  surface  S. 
N  is  the  source  term  that  arises  from  cylinderical  formulations. 


The  stale  vector  U  contains  five  or  more  variables  depending  upon  how  many  trace 
chemical  species  are  being  calculated  and  which  turbulence  model  is  used.  When  the 
k  e  turbulence  model  and  n  chemical  species  are  being  calculated,  U  is  defined  as: 


/  P  \ 


U  = 


pu 

pv 

pw 

E 

pk 

pt 

Pf\ 


\  Pfn  ) 
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The  flux  vector  P  is  defined  as 


p  =  Frl  +  Frj 


where  i  and  j  are  unit  vectors  pointing  in  the  axial  (x)  and  the  radial  (r)  directions, 
respectively.  Fx  and  Fr  are  the  axial  and  radial  components  of  the  flux  vector  defined 
below. 


Ft  = 


/  pu  \ 

puu  +  P  +  Txx 
puv  -f  Txr 
puw  -f-  Txg 

u(E  -f  P)  +  UTXX  +  VTrx  +  WTgx  -f  qx 
pku  —  Tk(dk/dx) 
peu  —  Tt(de/dx) 
pfiu  -  Th(dfi/dx) 

phu  -  rh(df2/dx) 

pfnu  -  Tu(dfn/dx) 


Fr 


f  pv  > 

pvu  +  Trx 
pW  +  Trr  +  P 

pVW  +  Trg 

V{E  +  P)  +  UTxr  +  VTrr  +  WTffr  +  qr 
pkv  —  rk{dk/dr) 
ptv  —  Y({dt/dr ) 
pfiv-rh(df,/dr) 
pf2v  -  T h(df2/dr) 


\  Pfnv-rjn(dfn/dr)  ) 
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The  source  term  N  arises  in  the  cylinderical  formulation  of  the  equations  and  is 
defined  below. 


N  = 


/  ° 

0 

( — pww  —  P  —  T6Q)/r 
(pvw  +  Trg)/r 
0 

pa  -  Gk 

~{Ct\  tGk  —  Ct2  p  c2)/ k 
0 
0 


V  0  ) 


The  definitions  of  some  of  the  terms  and  symbols  used  above  are  defined  as  follows: 

u  =  x- component  of  velocity 
v  =  r- component  of  velocity 
w  =  swirl-component  of  velocity 
e  =  specific  internal  energy 
p  —  mass  density 
P  =  pressure 
k  =  turbulent  kinetic  energy 
t  =  turbulent  energy  dissipation  rate 


E  =  pe  +  ~(u2  -f  v2  +  w2) 

it 
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— * 
V 

9* 

9r 


1  drv  du 
r  dr  +  dx 

-*  -y  «4 

=  UZ  +  VJ  + 

=  -r  E(dT/dx) 
=  -rE(dT/dr) 


I’m  is  the  effective  absolute  viscosity  defined  as  the  sum  of  the  laminar  and  turbulent, 
viscosities. 

Tm  =  f*i  + 

VE  is  the  effective  thermal  conductivity  defined  as  the  sum  of  the  laminar  and  tur¬ 
bulent  viscosities. 

T  e  =  fc(  +  fct 


F jt  is  the  effective  diffusion  coefficient  for  species  i  defined  as 


F  /,  =  Tm/o-/, 


where  cry,  is  the  schmidt  number  for  chemical  species  i. 

There  are  two  options  for  turbulence  models:  the  Baldwin- Lomax  model  and  the 
k-e  two-equation  model.  In  the  Baldwin  Lomax  model  the  turbulent  eddy  viscosities 
are  calculated  from  algebraic  relations  based  upon  local  quantities  of  the  flow  field. 
When  the  k-e  turbulence  model  is  used,  two  transport  equations  in  Equations  1  are 
solved.  The  pertinent  constants  and  relations  for  this  model  are  as  follows*. 

and  the  constants  are  defined  as 

Vt  = 

r*  = 

r,  = 

Cn  = 

C, a  = 
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C M  =  0.09 

(Tk  =  1.0 

<re  =  1.3 


The  values  of  Cc\ ,  C(2,  <?k,  &t,  and  may  be  specified  by  the  user.  The  above  values 
are  the  defaults. 

Equation  1  represents  a  very  simple  idea — the  time  rate  of  change  of  mass,  momen¬ 
tum,  and  energy  within  an  arbitrarily  chosen  volume,  V,  is  equal  to  the  apparent  flux 
of  these  quantities  inward  through  the  surface,  5,  surrounding  the  volume  plus  the 
production  of  these  quantities  within  the  volume.  The  finite  volume  method  consists 
of  breaking  the  flow  field  up  into  a  large  number  of  quadrilateral  finite  volume  cells 
to  which  the  integral  equations  are  applied. 


2.2  Thermodynamics 


The  equations  in  the  previous  section  must  be  supplemented  by  an  equation  of  state 
to  calculate  P  and  T  from  p  and  e.  The  DUMPSTER  code  offers  three  options  for 
equations  of  state:  1)  perfect  gas,  2)  thermally  perfect  gas,  3)  and  Tannehill’s[24] 
curve  fits  for  equilibrium  air. 


In  addition  to  calculating  P  and  T,  it  is  necessary  to  calculate  several  other  thermo¬ 
dynamic  variables  which  are  needed  to  calculate  the  numerical  fluxes.  The  first  two 
are  7  and  7.  For  perfect  gas  flows  7  and  7  are  identical  and  are  equal  to  the  ratio 
of  specific  heats,  For  real  gases,  however,  the  ratio  of  specific  heats  no  longer 
has  any  physical  significance  and  the  quantities  7  and  7  are  defined  to  simplify  the 
calculation  of  numerical  fluxes.  In  particular  we  use  Vinokur’s  definition [25], 


where  c,  is  the  isentropic  speed  of  sound.  For  numerical  fluxes  based  on  Roe’s  flux 
function  Vinokur  has  defined  the  alternative  variables  k  and  \  which  are  the  partial 
derivatives  of  pressure  P  with  respect  to  internal  energy  per  unit  volume  e  =  pe  and 
density  p. 
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K 


X 


dP 


de 


p 


dP_ 

dp 


i 


2.2.1  Perfect  Gas 


The  perfect  gas  equation  of  state  is  expressed  with  the  following  equations: 


P  =  (7  —  1 ) pe  =  *yRT 
T  =  e/Cv 


The  equation  of  state  can  be  altered  for  different  gases  by  specifying  the  ratio  of  spe¬ 
cific  heats,  7,  and  the  gas  constant,  R.  The  gas  constant  for  a  gas  of  a  given  molecular 
weight,  M,  is  calculated  from  the  universal  gas  constant  R  =  8.31434J/mole-K  using 
the  formula  R  —  R/M.  The  five  thermodynamic  variables  are  given  by  the  following 
formulas. 


7  =  constant 

7-1 

0 

R 

7-1 


2.2.2  Thermally  Perfect  Gas 


The  second  option  for  the  thermodynamic  equation  of  state  is  a  thermally  perfect 
gas.  For  this  gas  the  specific  internal  energy  (e)  is  not  linearly  proportional  to  the 
temperature.  In  general  the  internal  energy  and  temperature  relate  as  follows: 

e  =  e(T)  =  CV(T)T 
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where  CV(T)  is  the  temperature-dependent  specific  heat  capacity  at  constant  volume. 
For  a  perfect  gas  Cv  is  a  constant  and  the  temperature  can  be  easily  calculated 
from  the  internal  energy.  But  for  a  thermally  perfect  gas  where  Cv  is  a  function  of 
temperature,  calculating  temperature  from  internal  energy  requires  iteration. 

In  DUMPSTER  the  user  inputs  a  curve  fit  for  the  specific  heat  at  constant  volume 
as  a  function  of  absolute  temperature.  This  curve  fit  has  the  following  form: 

CV(T)  =  Clr2  -I-  c3T  +  c4 


To  avoid  an  iterative  process  every  time  temperature  is  decoded  from  internal  energy, 
a  mean  specific  heat  capacity  is  used  in  DUMPSTER  which  defined  as  the  ratio  of 
the  temperature  to  internal  energy. 

<7.  =  e/T  =  i  £  CV(T)  iT 


where  Ta  is  a  reference  temperature  arbitrarily  set  10  absolute  zero  ( T0  =  0).  An 
expression  for  the  mean  specific  heat  capacity  at  constant  volume  can  be  obtained  by 
integrating  the  curve  fit  function  for  specific  heat  capacity  with  respect  to  temperature 
as  follows: 


CV(T)  =  y)  TCl+l  +  yT2  +  c4T 


Temperature  is  decoded  by  dividing  the  internal  energy  by  the  mean  specific  heat. 


Cv  is  then  updated  at  the  new  temperature  for  use  in  the  next  time  step.  As  the 
solution  converges  so  does  the  value  for  Cv.  The  speed  of  sound  (c,)  used  for  the 
calculation  of  7  is  the  isentropic  speed  of  sound  defined  by 


Using  this  and  previously  given  definitions,  the  resulting  expressions  for  a  7  and  7 
for  a  thermally  perfect  gas  are 


7 

7 


,  R 
1  +  ~c 


,  R 
+  Cv 
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The  pressure  derivatives,  «  and  \i  are  needed  for  the  Roe  and  Harten-Yee  flux  func¬ 
tions  are  defined  as  follows: 


K  = 


X  = 


2.2.3  Equilibrium  Air  Curve  Fits 


The  third  option  for  the  equation  of  state  is  Tannehill’s  curve  fits  for  real  air.  With 
this  option  the  code  calls  a  subroutine  given  in  reference[24]  to  calculate  P  and  T 
from  p  and  e  for  each  finite-volume  cell  and  each  time  step.  This  routine  returns 
values  for  P ,  T,  and  c,  as  a  function  of  p  and  e  for  air  in  thermochemical  equilibrium 
up  to  25,000  K. 

Since  e0  =  0,  three  variables  can  be  calculated  directly  from  their  definitions  as 
follows: 


7  = 


7  = 


Cv 


pe 


e 


T 


The  pressure  derivatives,  k  and  must  be  calculated  using  data  numerical  differen¬ 
tiation  using  data  from  the  real  gas  curve  fit  routines.  To  evaluate  k  and  start 
with  the  chain  rule 


8P  = 


8p  =  k  8e  +  \8p 


and  the  expression  for  speed  of  sound  in  terms  of  k  and  \  as  described  by  Vinokur[25, 
Eq.  5]  (where  h  is  the  enthalpy) 


c]  =  X  +  *h 
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Eliminating  x  from  the  two  equations  above  and  solving  for  k  yields 

„  =  SP  -  c]6p 
Se.  —  hSp 

_  +  Sp,e  +  Se)P(p,  e)  -  c\Sp 

Se  —  hSp 


This  expression  may  be  used  to  calculate  k  from  any  combination  of  Sp  and  Se.  After 
some  experimentation,  Se  =  2h6p  wa s  chosen  since  it  seems  to  avoid  most  spurious 
results  near  discontinuities  in  the  curve  fits.  The  other  pressure  derivative,  x,  is  then 
calculated  from  x  =  —  tth. 


2.3  Transport  Properties 


The  effective  viscosity  I'm  and  the  effective  conductivity  T#  are  defined  as  the  sum 
of  the  laminar  and  turbulent  components  as  follows: 

I'm  =  Pi  +  Pt 
Te  =  ki  +  kt 


The  diffusion  coefficients  are  computer  from  T/v*  using  Schmidt  numbers.  These 
quantities  are  calculated  within  DUMPSTER  using  various  user  specified  methods  as 
described  below. 


2.3.1  Viscosity 


The  laminar  viscosity  is  calculated  using  one  of  three  methods.  The  first  method  is 
to  set  the  laminar  viscosity  to  a  constant.  The  second  method  is  to  assume  that  the 
viscosity  depends  only  on  temperature  and  that  the  dependence  may  be  approximated 
by  the  following  formula. 

cxTc 2  +  c3T  +  c4 

Pi  =  - - 

+  C6 

Sutherland’s  relation  for  viscosity  is  obtained  by  setting  the  constants  in  the  above 
formula  as  follows  where  the  temperature  ( T )  is  in  units  of  degrees  Kelvin  and  the 
viscosity  is  in  units  of  Kg/{m  ■  s). 

c,  =  1.451  x  1(T6 
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c2  =  1.5 
c3  =  0 
c4  =  0 
c5  =  1 
Cg  =  110 


The  laminar  viscosity  may  also  be  calculated  using  the  equilibrium  air  curve  fits  of 
Tannehill[26|.  A  viscosity  is  returned  when  the  Tannehill  routines  are  called  with  the 
density  and  internal  energy. 


2.3.2  Conductivity 


There  are  three  methods  for  calculating  the  thermal  conductivity.  One  method  is  to 
calculate  the  thermal  conductivities  from  the  viscosities  using  Prandtl  numbers  as 
follows: 


k,  = 


HiCp 

Pri 


kt  = 


f*tCp 

Prt 


where  Prt  and  Prt  are  the  laminar  and  turbulent  Prandtl  numbers  and  Cv  is  the  heat 
capacity  at  constant  pressure. 


The  thermal  conductivity  may  also  be  specified  as  a  polynomial  function  as  shown 
below. 

ciT02  +  c3T  +  C) 

ki  =  ■■  -  __ 

C5T  +  Cg 


The  laminar  conductivity  cooresponding  to  Sutherland’s  viscosity  for  air  and  a  Prandtl 
number  of  0.71  is  obtained  by  setting  the  above  constants  as  follows: 


=  2.053  x  10~3 
=  1.5 
=  0 


Cl 

C2 

C3 
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c4  =  0 
c5  =  1 

=  110 


The  third  method  for  calculating  the  thermal  conductivity  is  to  use  the  equilibrium 
air  curve  fits  of  Tannehill.  The  turbulent  conductivity  is  always  calculated  using  the 
turbulent  Prandtl  number  and  the  turbulent  viscosity. 


2.4  Boundary  Conditions 


The  boundary  conditions  are  an  important  aspect  of  the  problem  formulation.  At 
each  boundary  of  the  flow  domain,  boundary  conditions  must  be  applied  based  on 
the  physical  nature  of  the  boundary.  The  DUMPSTER  code  has  options  for  several 
types  of  boundary  conditions  as  described  in  the  following  subsections. 


2.4.1  Free-Slip  Walls 


When  an  in  viscid  analysis  is  being  performed,  or  when  a  known  stream  surface  is 
chosen  as  a  boundary  of  the  computed  flow  domain,  a  free-slip  boundary  condition  is 
specified.  The  boundary  condition  is  simply  that  there  is  no  flow  through  the  wall. 
This  means  that  the  component  of  velocity  normal  to  the  surface  is  zero  and  that 
the  component  of  velocity  tangent  to  the  surface  is  allowed  to  change  as  required 
by  the  governing  equation.  The  presence  of  a  free  slip  wall  also  affects  the  pressure 
and  temperature  fields.  If  the  wall  is  curved,  a  non-zero  pressure  gradient  is  needed 
normal  to  the  wall  to  turn  the  flow  in  the  direction  that  the  wall  is  curving.  This 
pressure  gradient  may  be  calculated  from  the  normal  momentum  equation  at  the 
wall  and  used  to  set  the  wall  pressure.  However,  in  the  analysis  of  dump  combustors 
and  nozzles,  the  viscous  nature  of  the  flow  near  actual  walls  is  generally  so  significant 
that  a  no-slip  boundary  condition  is  necessary.  The  free-slip  boundary  condition  then 
used  only  along  stream  surfaces,  such  as  planes  of  symmetry,  which  are  not  curved. 
The  gradients  of  pressure  and  temperature  normal  to  free-slip  walls  are,  therefore, 
assumed  to  be  zero. 
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2.4.2  No-Slip  Walls 


At  no-slip  walls  the  velocity  in  the  boundary  cells  is  set  to  the  negative  of  the  velocity 
components  in  the  corresponding  internal  cell  such  that  the  average  velocity  is  zero 
at  the  wall.  The  pressure  gradient  is  assumed  to  be  zero.  The  temperurature  in  the 
boundary  cell  is  set  so  that  the  specified  wall  temperature  results  when  the  internal 
cell  and  the  boundary  cell  temperatures  are  averaged. 


2.4.3  k-t  Turbulence  Model  Boundary  Conditions 

The  k-t  turbulence  equations  require  values  for  k  and  t  at  the  boundaries.  On  inflow, 
outflow,  symmetry,  and  free-slip  boundaries  these  quantities  are  set  in  the  boundary 
cells  assumming  zero  gradients. 

At  no-slip  boundaries  the  wall  function  is  used  to  calculate  the  wall  shear  stress 
and  calculate  the  values  for  k  and  t  at  the  center  of  the  internal  cell  next  to  the 
no-slip  wall.  Resolution  of  turbulent  boundary  layers  requires  a  very  finely  spaced 
computational  mesh  near  the  walls.  To  partially  alleviate  the  requirements  for  fine 
mesh  the  flow  near  the  wall  may  be  approximated  analytically  using  wall  functions  in 
DUMPSTER.  Wall  functions[22]  are  used  to  obtain  the  wall  shear  stresses.  As  shown 
below  a  logarithmic  function  relates  the  wall  shear  stress  to  the  velocity  at  the  center 
of  the  internal  cell  next  to  the  boundary. 


In  the  application  of  wall  functions  it  is  assumed  that  a  constant  shear  stress  region 
exists  in  the  boundary  layer  in  which  the  compressible  law-of-the-wall[27]  given  below 
is  valid. 

2l  =  I,n|W)+5.l 


The  shear  velocity  uT  is  defined  as 


and  u.  is  defined  as 


“•  =  fe)arcsi"{^©} 


<!>  =  hf-Ml  (l  +  Ml 
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In  the  above  equations  k  is  the  von  Karman  constant  with  a  value  used  here  of  0.4], 
Me  and  ue  are  the  Mach  number  and  velocity  tangent  to  the  wall  at  the  edge  of  the 
boundary  layer.  pu  and  /q  are  the  density  and  Iaininar  viscosity  of  the  fluid  at  the 
temperature  of  the  wall,  r/  is  a  length  scale  equal  to  the  distance  from  the  wall  to  the 
center  of  the  internal  cell  next  to  the  boundary  where  u  is  defined. 


2.4.4  Supersonic  Inflow 


The  supersonic  inflow  boundary  condition  requires  specification  of  all  variables.  This 
boundary  condition  is  only  appropriate  where  the  Mach  number,  based  on  the  velocity 
normal  to  the  boundary,  is  greater  than  one. 


2.4.5  Subsonic  Inflow 


At  a  subsonic  inflow  boundary  the  total  pressure  (Pj-),  total  temperature  (Tj),  and 
flow  angles  are  specified  as  a  function  of  the  radius  r.  Along  the  inflow  boundary 
the  total  velocity  magnitude  (|V|)  is  extrapolated  from  within  the  flow  field.  Using 
the  isentropic  relation  below  the  local  Mach  number  (M)  is  calculated  for  the  given 
velocity  magnitude  and  total  temperature. 


M  = 


2  q 


7(7  -  1)CvTt  -  (7  -  1)? 


where  <7  is  defined  as  q  =  \(u2  +  v2  +  w2).  Using  the  isentropic  relations  below,  static 
pressure  and  temperature  can  then  be  calculated.  Prescribed  flow  angles  (a,)  are  used 
to  set  the  boundary  velocity  components  given  the  extrapolated  velocity  magnitude. 

P  =  pT(\ +^{7-  \)M2y^ 

T  =  7V(l+^(7-l)M2)  ‘ 

U  =  QX|U| 

V  =  C*r|U| 

w  =  a*|U| 
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2.4.6  Supersonic  Outflow 


Along  a  supersonic  outflow  boundaiy,  the  flow  variables  in  the  boundary  cells  are  com¬ 
pletely  determined  by  the  upstream  flow  conditions.  Strictly  speaking,  this  boundary 
condition  is  only  appropriate  when  the  Mach  number,  based  on  the  velocity  normal  to 
the  boundary,  is  greater  than  one.  Its  use  is  acceptable  for  boundary  layer  flow  when 
the  boundary  layer  is  thin  and  the  flow  outside  of  the  boundary  layer  is  supersonic. 


2.4.7  Subsonic  Outflow 

The  st  atic  pressure  on  the  outflow  boundary  is  specified  at  the  upper  nozzle/combustor 
wall.  The  actual  outflow  boundary  static  pressure  must  vary  in  the  radial  direction 
when  axial  swirl  is  present.  According  to  the  simplified  radial  momentum  equation: 

dP  1  2 

"ST  =  «PW 


This  equation  is  integrated  to  obtain  the  outflow  boundary  static  pressure  distri¬ 
bution.  The  values  for  pk ,  pe,  and  pf,  a*  *he  outflow  boundary  are  obtained  by 
zero-order  extrapolation  from  the  interior. 


2.4.8  Interzone  Boundary  Condition 


At  the  boundary  between  adjacent  zones,  the  conservative  variable  data  are  simply 
transferred  from  the  internal  cells  of  one  zone  to  the  boundary  cells  of  the  adjacent 
zone. 


2.5  Initial  Conditions 


In  DUMPSTER  there  are  three  options  for  specifying  initial  flow  conditions  as  de¬ 
scribed  below.  For  each  computational  zone  one  of  the  following  methods  may  be 
used. 

•  Uniform  Conditions  The  flow  field  is  set  to  uniform  conditions  as  specified 
for  this  zone. 
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•  One-Dimensional  Solution  The  total  pressure  and  total  temperature  at  the 
inflow,  the  area  variation,  and  the  effective  throat  area  are  specified.  From  this, 
a  one-dimensional  isentopic  solution  for  u,  P,  T  and  p  is  calculated.  The  cross- 
flow  velocity  component  v  is  set  so  that  the  flow  is  parallel  to  the  walls.  The 
swirl  component  of  velocity  is  uniform. 

•  Custom  The  user  may  provide  FORTRAN77  coding  that  generates  custom 
initial  conditions. 
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3  DISCRETIZATION 


The  computational  domain  is  divided  into  one  or  more  computational  zones.  As 
shown  in  Figure  2  these  zones  are  stacked  in  the  r  (radial)  direction.  Multiple  zones 
allow  easy  mesh  generation  for  a  wide  range  of  configurations  of  dump  combustors 
and  nozzles.  The  domain  within  each  zone  is  discretized  into  a  grid  of  finite-volume 
cells  as  shown  in  the  figure.  Within  each  zone  the  cells  are  indexed  by  two  parameters: 
i  and  j. 

An  individual  finite-volume  cell,  with  indices  i  and  j ,  is  shown  in  f  igure  3.  Applying 
the  Navier-Stokes  equations  in  integral  form  to  this  volume  gives 

=  -(A  P-S+Dj  P-S)  +  Nij 


where  A,j  is  the  mean  value  of  U  in  cell  i,j  and  DiPS,  for  example,  represents  the 
difference  of  the  fluxes  through  opposing  i-faces  of  the  cell. 

The  time  derivative  is  approximated  using  backward-in-time  differencing  where  n  is 
the  old  time  level  and  n  +  1  is  the  new  time  level. 


=  -(Di  P  S  +  Dj  P-S)n+'  +  /V"+1 


C2) 


where 


Figure  2:  Multiple  Computational  Mesh  Zones. 
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®  1*1/2  J 


i 


Figure  3:  Finite- Volume  Cell. 

For  the  approximation  of  the  flux  through  a  surface  the  inviscid  and  diffusion  terms 
of  the  flux  vector  arc  considered  separately. 

P-s  =  P-sinv-  +  Psd'fi- 


These  terms  are  then  evaluated  in  a  manner  consistent  with  the  predominant  nature 
of  the  equations  in  the  limit  as  Re  — »  oo  (hyperbolic)  and  Re  —*  0  (parabolic)  (i.e., 
upwind  differencing  for  the  inviscid  terms  and  central  differencing  for  the  viscous 
stress  and  heat  flux  terms.  Three  flux  functions  are  have  been  implemented  into 
DUMPSTER: 


•  Steger  and  Warming’s  flux-vector-splitting[28], 

•  Harten  and  Yee’s[29]  second  order  TVD2  scheme,  and 

•  Roe’s  Approximate  Reimann  solver[30]. 


Each  of  the  above  flux  functions  is  described  in  the  following  subsections.  The  diffu¬ 
sion  terms  are  evaluated  using  standard  central  differences [31]. 

2Total  Variation  Diminishing 


21 


Figure  4:  Mathematical  Characteristics  of  One-Dimensional  Subsonic  Flow. 


3.1  Steger  and  Warming  Flux  Function 


Flux-vector  splitting  was  developed  by  Steger  and  Warming[28]  and  Van  Leer[33]  as 
a  way  to  upwind  difference  the  Euler  equations  in  regions  of  subsonic  flow.  If  the  flow 
normal  to  the  surface  of  a  cell  is  supersonic,  the  theory  of  characteristics  tells  us  that 
the  flow  field  at  the  surface  depends  only  on  the  solution  upstream  of  the  surface. 
In  this  case,  the  flux  at  the  surface  may  simply  be  evaluated  using  the  solution  in 
the  cell  upstream  of  the  surface.  The  numerical  flux  ( F )  at  the  i  +  l/2,y-face  of  the 
i,  j-cell  is  given  as 


P  f  F(ttj)  if  u' >c. 

,+1/2j~\  F(Ui+u)  if  u'<c. 


(3) 


Where  u'  is  the  velocity  component  normal  to  the  cell  face.  If  the  flow  is  subsonic, 
however,  the  flux  through  the  t  +  l/2,  j  surface  depends  on  the  solution  both  upstream 
and  downstream  of  the  surface.  Then  the  simple  upwinding  applicable  to  supersonic 
flows,  Equation  3,  is  no  longer  appropriate.  The  term  upwind  differencing  must 
now  be  defined  to  mean  differencing  in  a  manner  appropriate  with  the  mathematical 
characteristics  of  the  flow.  As  shown  in  Figure  4,  there  are  characteristics  running 
both  to  the  right  and  left  in  subsonic  flow  and  upwind  differencing  must,  therefore, 
use  backward  differencing  on  some  of  the  flux  and  forward  differencing  on  the  rest  of 
the  flux. 
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(4) 


In  flux- vector-splitting  methods,  the  flux  function  is  written  as 

F(U)  =  F+{U)  +  F~(U) 

so  that  the  flux-split  Jacobians,  A±  =  have  eigenvalues  satisfying 

Am  (A+)  >  0  for  all  m 
Am  (J4_)  <  0  for  all  m  (5) 

Flux  splitting  has  divided  a  flux  which  cannot  in  general  be  upwinded  into  two  fluxes 
which  separately  can  be  upwinded.  With  first-order  upwinding  the  flux  through  the 
i  -f  1/2,  j  surface  is 

Fi+tm  =  F*  (Vij)  +  F-  (W+u) 

Higher-order  accurate  methods  are  obtained  by  combined  extrapolation  and  interpo¬ 
lation  of  the  solution  to  the  surface  from  nearby  cells. 

v-  =  uw+|[a-«)(«j-t/i-.j)+(i+u)(wi+w-i'w)]  w 

U+  =  tfi+ij  +  !  1(1  -  ■*>)  (l/i+1  j  -  UiJrU)  +  (1  +  u)  (U.j  -  £/,+w)|  (7) 

Here  (f>  varies  the  differencing  between  first-order,  <f>  —  0,  and  second  order,  <f>  =  1. 
The  second  parameter,  u,  varies  the  differencing  between  fully  upwind,  u>  =  —  1,  and 
central  differencing,  u  =  1.  Using  MUSCL-differencing[32]  the  second-order  flux  is 
then 

£+./*,  =  F *  (tr)  +  F-  (u+) 

Equations  4  and  5  do  not  uniquely  define  the  split  fluxes  F±  (U).  In  fact,  there  are 
many  flux-splittings  that  have  been  proposed[28,  33,  34].  The  two  most  common 
flux-splittings  are  those  of  Steger  and  Warming[28]  and  Van  Leer[33].  Both  splittings 
were  originally  developed  for  ideal  gases  and  were  later  extended  for  use  with  real 
gases  [25,  12]. 

Steger  and  Warming  developed  their  splitting  based  on  the  homogeneity  property, 
F  =  AU,  of  the  flux  vector  with  a  thermally  perfect  gas.  They  wrote  the  flux 
Jacobian  as 

A  =  A+  +  A~ 

where 

A±  =  X~l  A±X 
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and  A*  is  the  diagonal  matrix  with  elements 

A*  =i(A„±|Am|) 

Here  the  matrices  A,  X,  and  X~l  are  the  diagonal  matrix  of  eigenvalues,  matrix  of 
eigenvectors,  and  its  inverse,  for  the  flux  Jacobian  A.  Therefore,  A+  and  A~  arc 
equal  to  the  flux  Jacobian  matrix  with  the  negative  and  positive  eigenvalues  set  to 
zero,  respectively.  The  split  fluxes  are  F ±  =  A*U.  These  matrix  multiplies  may  be 
performed  numerically  for  each  face  on  each  time  step. 

Significant  savings  in  computer  time  is  possible  by  eliminating  the  matrix  multiplies 
as  follows.  Consider  the  contribution  to  the  total  flux  from  each  of  the  following  three 
unique  eigenvalues  separately. 

A,  =  u' 

A2  =  u'  +  c, 

A3  =  u  -  c. 


The  Aj  eigenvalue  is  repeated  so  that,  for  flows  without  species  transport, 

A  =  diag  { A , ,  A2,  Ai,  Aj,  A3} 

The  separate  contributions  from  each  eigenvalue  are  obtained  by  defining  diagonal 
matrices  with  all  but  the  desired  eigenvalue  set  to  zero. 

Aj  =  diag  {Aj,0,  AI,  A,,0} 

A2  =  diag  {0,  A2, 0,0,0} 

A3  =  diag  {0,0,0, 0,  A3} 

The  fluxes  are  then  Fm  =  X-1  AmX U .  Performing  the  above  operations  for  thermally 
perfect  gas: 

1 
u 
v 
w 

i(u2  +  u2  +  u,>)  +  pe-^ 
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•% 


F3 


A ip 
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w 
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1 

u  +  c,nx 
v  +  c3nr 
w  +  c„n2 
H  +  c,u' 

1 

u  -  csnx 
v  -  c„nT 
w  —  c,n2 
H  -  csu' 


(9) 


(10) 


u'  is  defined  as  u'  =  nxu  +  nT v  +  nzw.  Vinokur[25]  used  an  alternative  approach  to 
derive  the  Steger  and  Warming  fluxes  and  found  that  the  above  formulas  also  apply 
to  real  gases,  provided  that  7  is  defined  as 


The  split  fluxes,  F±,  are  obtained  from  the  Fm  by  splitting  according  to  the  sign  of 
the  eigenvalue.  For  —  ca  <  u'  <  0  we  therefore  have 

F+  —  F2  and  F~  —  Fj  +  F3. 

These  fluxes  are  relatively  inexpensive  to  calculate  compared  with  the  formulation 
requiring  numerical  matrix  multiplies. 

In  general,  some  form  of  flux  limiting  is  required  to  avoid  oscillations  at  shock  waves 
and  where  large  second  derivatives  of  the  flow  field  variables  exist.  Limiting  is  im¬ 
plemented  by  locally  varying  the  difference  quotients  in  the  MUSCL-differencing  in¬ 
terpolation  formulae  given  above  in  relation  to  the  local  gradients  of  the  solution. 
If  properly  implemented  a  limiter  will  yield  sharp  monotonic  shock  waves.  In  this 
investigation,  the  limiter  is  based  on  the  second  difference  of  pressure.  In  particular 

4>i+ 1/2  =  MAX  ^1  —  Cum-  ,  O.Oj 

where  AAp  is  normal  to  the  i  +  1/2,  j  cell  face  (in  this  case,  the  t-direction).  The 
differencing  parameter,  a;,  remains  unchanged.  This  limiting  reduces  the  form  of 
differencing  toward  first-order  upwind  differencing  in  regions  of  strong  pressure  gra¬ 
dients  (such  as  shock  waves).  The  spacial  order  of  accuracy  remains  second-order 
everywhere  unless  <f>  is  actually  zero. 
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3.2  Harten  and  Yee  Flux  Function 


In  the  eaxly  1980’s  Yee,  Warming,  and  Harten[29],  and  others[35,  36,  37]  presented 
upwind  biased  schemes  for  nonlinear  scalar  equations  for  which  the  total  variation  of 
the  solution  always  diminished  as  the  solution  proceeds. 

TV  (f/n+t)  <  TV  (Un) 


where 


TV(V)  =  £  |l/,+1  -  V. | 


Schemes  which  are  total  variation  diminishing  (TVD)  are  guaranteed  not  to  generate 
spurious  oscillations  and  are  therefore  robust  for  applications  involving  shock  waves. 
When  the  same  scheme  is  applied  to  nonlinear  systems  it  is  no  longer  TVD,  but  it  may 
still  exhibit  good  accuracy  and  robustness  for  problems  involving  strong  shock  waves. 
A  finite-volume  version  of  Harten  and  Yee’s  scheme  is  included  in  the  DUMPSTER 
code.  The  scheme  is  obtained  by  using  the  Harten-Yee  flux  function  described  in  this 
section. 

The  Harten-Yee  flux  function  is  written  as  a  central  difference  flux  plus  a  dissipation 
term. 

^.+  1/2, j  =  2  [F.+1.J  +  F'i.j  +  ('^4^>|)1  +  i/2,j]  (11) 

The  dissipation  operators  are  defined  as  follows: 

(^4)i+i/2  —  9\  +  9t+\  ~  'l'  (‘Y+1/2  +  1  '+1/2)  <*1+1/2  (12) 

where 

ai+l/2  =  (Xa),+\/2  Wi+l.j  ~  Ui,j)  (13) 

and,  in  the  standard,  form  'k(m(2)  =  l^ml-  A  modified  form  o<"  ^(r)  will  be  discussed 
latter  in  this  section. 

A  key  to  the  definition  of  the  Harten-Yee  flux  is  the  eigensystem  of  the  flux  Jacobian 
matrix,  A  =  The  components  of  the  eigensystem  are  the  diagonal  matrix  of 
eigenvalues  A,  the  matrix  of  right  eigenvectors  X,  and  its  inverse  (A")-1.  In  the  above 
equations,  4>  is  the  dissipation  term  which  will  reduce  the  scheme  to  Roe’s  first-order 
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upwind  scheme  in  regions  of  steep  gradients.  The  definition  of  the  flux  is  completed 
by  the  relations 


gt  =  S  max  [o,min  (ff,>i/2|a«+i/2|,5<7i-i/2ai-  1/2)] 
S  =  sign(ai+i/2) 

<*»+ 1/2  =  2^  (^*+1/2) 


and 


T.+1/2 


■! 


g'  +  l~gt 
0 


if  <*i+1/2  ±  0 

if  q«+i/2  =  0 


To  avoid  excessive  roundoff  error,  Tj+1/2  in  the  above  equation  is  actually  set  to  zero 
any  time  a,+  i/2  is  near  zero  (  — 10~22  <  c*,+i/2  <  10-22). 

This  scheme  is  theoretically  second-order  accurate  in  space.  In  regions  where  the 
solution  is  smooth  the  dissipation  term  is  turned  off  and  the  scheme  reduces  to  a 
central  difference  scheme.  Conversely,  for  regions  with  steep  gradients,  such  as  shock 
waves,  it  reduces  to  a  first-order  scheme  to  avoid  aphysical  oscillations  in  the  solution. 
The  first-order  scheme  to  which  it  reduces  is  the  flux-difference  splitting  of  Roe[30]. 

Since  the  Harten-Yee  flux  function  is  based  on  Roe’s  flux-difference  splitting,  it  shares 
certain  abnormalities  with  Roe’s  scheme.  In  particular,  Roe’s  scheme  does  not  nec¬ 
essarily  satisfy  the  entropy  inequality  and  may,  therefore,  yield  aphysical  solutions 
such  as  expansion  shock  waves[35,  36,  37],  To  overcome  this  weakness,  'F  is  modified 
as  follows: 

—  ~  "F  \J  z  mm  d”  cm^ 


Far  away  from  zero  each  element  of  'F  approaches  the  absolute  value  of  the  corre¬ 
sponding  element  of  z,  as  it  should.  When  zmm  approaches  zero,  however, 
approaches  the  positive  number  em. 

The  values  of  em  are  calculated  based  on  user  defined  constants  and  the  type  of 
characteristic  equation  which  corresponds  to  Am.  If  the  characteristic  transports 
entropy 

Cm  =  [C5  +  (M  +  C,)] 


If  the  characteristic  transports  acoustic  waves 

tm  =  [e*i  +  (|u|  +  c,)] 
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Finally,  if  the  characteristic  transports  momentum  components  tangent  to  the  surface 


t-m  —  [c7  +  (M  +  C,)] 


Different  constants  were  provided  for  each  type  of  characteristic.  A  flow  problem  with 
shocks  may  require  a  large  value  of  <%  to  avoid  an  expansion  shock,  but  a  small  value 
for  c-r  to  avoid  excessive  smearing  of  the  boundary  layer. 

The  CPU  time  for  the  matrix  multiplies  in  Equations  11  and  13  would  normally 
increase  quadratically  with  the  number  of  transport  equations.  In  DUMPSTER  how¬ 
ever,  the  matrix  multiplies  are  performed  analytically  whenever  possible,  and  the  CPU 
time  increases  nearly  linearly  with  the  number  of  added  specie  transport  equations 
and  k-t  equations. 

To  evaluate  the  variables  at  the  cell  surfaces,  arithmetic  averaging  of  the  cell  centered 
values  is  used.  This  kind  of  averaging  has  the  advantage  of  computational  simplicity 
and  can  be  easily  extended  for  problems  in  thermochemical  nonequilibrium.  Other 
kinds  of  averaging  can  also  be  used,  but  they  are  usually  more  complicated.  Arith¬ 
metic  averaging  takes  orm 

1  / 

A+i/2  =  2  pi+l ' 

The  same  procedure  is  used  to  calculated  u,  v,  iu,  and  a  at  the  surface. 


3.3  Roe’s  Flux  Function 


Roe’s  flux  function  is  obtained  when  the  ^’s  in  the  discussion  fo  the  Harten-Yee  flux 
functions  are  set  to  zero.  The  method  is  first-order  accurate. 
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4  SOLUTION  PROCEDURES 


4.1  Line  Gauss-Seidel  Method 


The  Line  Gauss-Seidel  (LGS)  method  is  an  unfactored  method  based  on  the  work 
by  Thomas  and  Walters[17],  MacCormack[16],  and  Chakravarthy[18].  This  type  of 
algorithm  has  proven  to  be  an  efficient  relaxation  procedure  for  steady-state  flow 
fields.  The  algorithm  used  in  this  investigation  is  an  improvement  on  the  algorithm 
presented  by  Peery  and  Imlay[38,  39].  A  brief  description  of  the  basic  algorithm  is 
given  below. 

For  conciseness  we  consider  only  the  calculation  of  the  flux  at  the  i  +  l/2,j  cell  face  as 
shown  in  Figure  3.  The  left  side  of  the  face  refers  to  the  side  that  of  the  i,  j  cell,  and 
the  right  side  of  the  face  refers  to  the  side  of  the  i  +  l,j  cell.  For  the  approximation 
of  the  flux  through  a  surface  the  inviscid  and  diffusion  terms  of  the  flux  vector  are 
considered  separately. 


P  S  =  p-s  'mv-  +  P  S  diff- 


The  total  flux  vector  is  a  function  of  the  solution  in  nearby  cells 

^  ■^1+1/2  =  Ui.Ji  £4+l,J)  ^t+2,j) 

The  flux  is  evaluated  at  the  new  (unknown)  time  level  and  linearized  as  follows: 


p  Q  n+1 
r  'D  i+\/2 


P'S  i+ 1/2  + 

>+l 


^1  +  1/2 ,j  T  A1+1/2,j  SUij  +  ^i+l/2,j  ^»+lj  T  A,+lyf2j  6Ui+2j 


(14) 


where 


l‘-l 

i+l/2 

i+l/2,j 

1 

6" 

1 

I* 

Sfi+1/2 

dU, 

,<  +  l 

dfow 

•+1/2.J 

dU,+\ 

l«+2 

O&ip 

dUi+ 2 
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The  equation  for  the  flux  at  the  Sj+ j/2  face  is 


Bi'Un  +  Br+1»  SUtj  +  6Uij»  +  Bjt'y,  (15) 


where 


Eji-i 

Aj+l/2  “ 

dUj.  i 

DJ  _ 

3/"+l/2 

a>j+ 1/2  “ 

at/, 

RJ  +  1  _ 

°ij+ 1/2  ~ 

dfj+ 1/2 

DJ+2  _ 

A.j+1/2  - 

5/JVl/2 

«^+2 

The  above  Jacobians  (A  and  B)  are  calculated  numerically  on  each  cycle.  Substituting 
Equations  14  and  15  into  the  discrete  conservation  equation  (Equation  2),  grouping 
terms  and  linearizing  the  source  term  (Nli})  yields  the  following  finite-volume  equa¬ 
tion. 


V.-J 


1  dNij 

-  4 - hi 

6t  dUij 


Di  K;|/2j  W-H  +  A  i+1/2j  +  A!*\,u  sv»u  +  A'*\,XJ  SUM,,\+ 

D,  \BijlwWj-,  +  +  Bj£,/2SUiJ+l  +  =  By  (16) 


where  the  residual  Rij  is  defined  as 

Rij  =  ~  (A  P  S  +  Dj  P  S  +  V,,j  Nij)  "  (17) 

The  residuals  approach  zero  as  the  solution  approaches  steady  state.  On  each  time 
step  the  solution  is  advanced  as  follows: 


1.  Calculate  residuals  based  upon  the  local  variation  of  the  solution  for  each  cell 
using  the  explicit  flux  balance  given  in  Equation  17. 
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Figure  5:  Linearized  Algebraic  System  of  Equations. 

2.  Calculate  the  change  in  the  solution  by  solving  the  block-linear  system  of  alge¬ 
braic  equations  shown  in  Figure  5.  The  rows  of  this  system  are  the  block-linear 
algebraic  relations  resulting  from  the  implicit  flux  balance,  Equation  16,  applied 
to  each  cell.  This  step  implicitly  calculates  the  change  in  the  solution  according 
to  global  variations  in  the  residual. 

2.  Update  the  solution  for  each  cell  using 

c;;1  =  c,"  +  (u  su,j 

where  J  is  a  local  relaxation  factor. 

The  first,  and  last  steps  of  the  above  algorithm  are  easy  to  implement.  The  second  step 
is  difficult,  since  it  involves  the  solution  of  block  banded  system  of  linear  algebraic 
equations,  as  shown  in  Figure  5.  The  coefficient  matrix  for  this  system  is  very  large, 
144,000x144,000  for  a  60x60  mesh,  but  is  sparse  and  well  structured.  Because  the 
system  is  so  large  it  is  impractical  to  solve  it  directly  and  most  investigators  resort  to 
approximately  factoring  this  matrix  into  simpler  matrices  such  as  block  tridiagonals. 
Unfortunately  the  error  associated  with  the  approximate  factorization  can  severely 
restrict  the  allowable  time  step,  and  hence  convergence  rate.  In  this  investigation  this 
system  was  solved  iteratively  using  a  modified  alternating  direction  line  Gauss-Seidel 
relaxation  method. 

When  the  linear  system  of  equations  is  solved  accurately  and  the  time  step  is  taken 
very  large,  this  procedure  approximates  a  Newton  method  for  the  system  of  nonlinear 
algebraic  equations  resulting  from  differencing  of  the  time  independent  Navier-Stokes 
equations.  As  with  any  Newton  method,  this  procedure  will  only  converge  if  the 
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chosen  initial  conditions  are  sufficiently  close  to  the  final  solution.  Other  similar 
methods[l7,  16]  have  avoided  these  convergence  problems  by  starting  with  small, 
nearly  explicit,  time  steps  and  increasing  the  time  step  as  the  solution  progresses. 
Unfortunately,  it  is  impossible  to  choose  a  priori  a  time  step  sequence  which  will 
yield  both  a  stable  solution  and  an  optimal  convergence  rate.  Generally,  when  a  new 
problem  is  attempted,  the  time  step  sequence  is  such  that  the  solution  diverges  or 
it  takes  many  times  the  optimal  number  of  cycles  to  converge.  As  a  result,  many 
attempted  runs  are  often  required  to  obtain  one  solution. 

These  problems  were  overcome  by  adapting  a  technique  for  global  convergence  used  in 
solving  other  systems  of  nonlinear  equations.  A  controlling  subroutine  automatically 
adjusts  the  relaxation  factor  the  time  step  size,  the  number  of  explicit  steps 

between  each  implicit  step,  the  relaxation  factor  for  the  Gauss-Seidel  iteration  of  the 
linear  system  The  technique  is  described  in  more  detail  by  Peery  and  Imlay[38], 

The  implicit  time  step  is  still  set  small  initially,  but  is  increased  exponentially  to  a  very 
large  value  ( CFL  «  106).  The  overhead  for  this  procedure  depends  on  the  problem, 
but  is  certainly  acceptable  considering  that  it  is  only  used  if  it  is  needed.  With  this 
technique  the  reliability  of  the  Navier-Stokes  solver  has  increased  dramatically  and 
converged  blunt  body  solutions  are  routinely  obtained  in  10  to  100  implicit  steps. 


4.2  Lower-Upper  Symmetric  Gauss-Seidel  Method  (LU-SGS) 


The  LGS  method  explained  above  proved  to  very  fast,  but  required  a  flux  function 
that  is  numerically  dissipative.  The  Lower-Upper  Symmetric  Gauss-Seidel  (LU-SGS) 
method  developed  by  Yoon  and  Jameson[20]  was  also  studied  in  this  work.  The 
LU-SGS  method  is  not  as  fast  as  the  LGS  method,  but  the  LU-SGS  method  is  more 
robust  and  reliable,  especially  when  the  flux  functions  have  low  numerical  dissipation. 
In  this  work  having  low  numerical  dissipation  is  important  for  accurately  simulating 
the  recirculation  regions  of  the  dump  combustor. 

The  TVD  flux  function  of  Yee  and  Harten  was  adopted  for  use  with  the  LU-SGS 
method,  since  it  exhibits  low  numerical  dissipation.  The  numerical  dissipation  inher¬ 
ent  with  flux  functions  can  overpower  the  natural  action  of  viscosity  in  shear  layers 
unless  very  fine  meshes  are  used.  This  is  especially  true  with  the  Steger  and  Warming 
flux  function.  The  Harten- Yee  TV1)  flux  function,  on  the  other  hand,  is  less  dissi¬ 
pative  and,  therefore,  requires  fewer  mesh  points  for  a  given  level  of  accuracy  than 
either  of  the  previous  flux  functions.  The  LGS  method  would  not  converge  with  the 
Harten  Yee  TVD  flux  function. 


The  LU-SGS  implicit  algorithm  requires  less  memory  than  the  LGS  algorithm.  With- 


out.  species  transport  or  the  k-t  turbulence  model  equations,  the  LGS  algorithm  re¬ 
quired  storage  for  187  variables  at  each  mesh  cell.  Under  the  same  conditions  the 
blJ-SGS  algorithm  requires  storage  of  27  variables  per  mesh  cell.  The  difference  is 
even  more  dramatic  when  the  k-t  turbulence  model  and  specie  transport  equations 
are  considered.  With  the  old  flux  functions,  the  added  speed  of  the  LGS  algorithm 
justified  its  inefficient  memory  use,  but  with  the  adoption  of  the  TVD  fluxes,  that  jus¬ 
tification  is  gone.  Another  advantage  of  the  LU-SGS  is  that  it  can  easily  be  extended 
to  chemically  reacting  flows  in  which  the  the  addition  of  chemical  species  increases 
the  computational  work  proportionally. 


The  flux  vector  is  a  function  of  the  solution  in  nearby  cells.  Consider  only  the  2  +  1/2,  j 
surface  of  the  i,  j  cell  as  shown  in  Figure  3. 


^•*^+1/2  —  f(Ui-l  J,  f/ij,  Ui+l,j,  Ul+2,j) 


(18) 


The  flux  is  evaluated  at  the  new  (unknown)  time  level  and  linearized  using  Yoon’s 
approximate  .Jacobians. 

(P-S  *  (P-5)?+1/2  +  +  4r+1/JJ«/,+u  (19) 

where 


1/2, j  — 

2  (^».j  r(Ai,j)) 

^i+l/2,j  = 

= 

dPij 

dUij 

and  r a  is  the  spectral  radius  (maximum  eigenvalue)  of  A.  Substituting  Equation  19 
and  the  corresponding  flux  for  the  /-direction  faces  into  the  discrete  conservation 
equation,  Equation  2,  yields  the  implicit  finite-volume  equation.  This  equation  may 
be  written  as  three  steps: 


1.  Calculate  the  residuals  for  each  cell  using  an  explicit  flux  balance. 

RtlJ  =  -(D,  P-S  +  D}  P-S)n  +  JV" 

The  residuals  will  approach  zero  as  the  solution  approaches  steady  state. 

2.  Calculate  the  change  in  the  solution  by  solving  the  block-linear  system  of  alge¬ 
braic  equations.  The  rows  of  this  system  are  the  block-linear  algebraic  relations 
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resulting  from  the  implicit  flux  balance  applied  to  each  cell. 


err  |  err  . 

+  W,  iJ 

Di  fiUij  "b  ^i+i/2,i  "b  Dj  [Btj+l/2  “b  Bi  J  +  j/2  SUij+i]  =  Ri,] 

Because  of  the  choice  of  Jacobians,  this  resulting  block  matrix  multiplying  6UtJ 
is  approximately  diagonal  when  the  source  Jacobian,  §77^,  is  zero. 


3.  Update  the  solution  using 


VS'  =  V?J  +  (rSUij 


where  £r  is  a  relaxation  factor.  DUMPSTER  starts  with  a  users  specified  £r 
near  one  and  updates  the  solution.  If  the  solution  at  time  level  n  +  1  results 
in  negative  pressures  or  temperatures  the  relaxation  factor  is  halfed  and  the; 
solution  at  n  -f  1  is  recalculated.  This  process  is  continued  until  all  pressures 
and  temperatures  are  non- zero  or  £r  reaches  a  specified  minimum. 


The  above  algorithm  consists  of  three  steps  which  must  be  done  on  each  time  step. 
The  first  step  is  an  explicit  step  which  calculates  the  residual  from  the  local  variation 
of  the  solution.  The  second  step  implicitly  calculates  the  change  in  the  solution 
according  to  global  variations  in  the  residual.  The  last  step  updates  the  solution. 

The  second  step  is  the  most  difficult  because  it  requires  the  solution  of  a  large  sparse? 
system  of  linear  equations.  The  LU-SGS  algorithm  approximately  solves  this  system 
using  two  sweeps  of  a  point  Gauss-Seidel  relaxation.  Assume  the  iteration  sweeps 
first  in  the  direction  of  increasing  i,j  indices  and  then  in  the  direction  of  decreasing 
j  indices.  The  increasing  i,j  iteration  is  performed  by  applying  the  equation 

iir+/3 ,rA+0ra),+{w)iJ 

(20) 


Wj  -  Atl/2,jSUi-hj 


=  flo¬ 


at  each  internal  cell.  The  decreasing  i,j  iteration  is  performed  by  applying  the  equa¬ 
tion 


+ /3rA  +  I  + 


(2*) 


This  scheme  may  be  written  as  an  approximate  lower-upper  (LU)  factorization  by 
subtracting  Equation  20  from  Equation  21  and  rearranging  to  remove  61  f*  from  the 


M 


left  hand  side. 


St 


(1  +  fir A  +  fir B)  I  + 


f-1  1 

vw/J 


V«.i 


ft 


(1  +  0rA  +  0rB)l  + 


(w)  «:* 

X  '  «.jJ 


(22) 


0  is  a  coefficient  used  to  artificially  increase  the  spectral  radius  and  thereby  increase 
the  numerical  stability.  The  coefficient  is  usually  set  to  1.0.  Increasing  it  above  1.0 
will  normally  reduce  the  convergence  rate.  The  equation  above  replaces  Equation  21. 
It  yields  the  same  results  as  Equation  21,  but  requires  fewer  floating  point  operations 
and  less  memory. 

The  source  Jacobian,  is  generally  not  diagonal.  However,  for  this  investigation  it 
was  assumed  to  be  zero  with  no  adverse  consequences.  With  this  assumption  the  LU- 
SGS  algorithm  requires  no  block  matrix  inversions  and  is,  therefore,  very  inexpenxive 
per  time  step  compared  to  the  LGS  algorithm.  This  advantage  is  counter  balanced 
by  the  comparatively  slow  convergence  rate  of  the  LU-SGS  algorithm  which  may  take 
thousands  of  time  steps  as  compared  to  hundreds  of  steps  for  the  LGS  algorithm. 
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5  COMPUTER  CODE 


DUMPSTER  is  written  in  the  FORTRAN77  language  and  may  be  compiled  on  almost 
any  computer.  It  runs  in  batch  mode,  reading  and  writing  to  several  input/output 
files.  It  consists  of  one  main  program  module,  57  major  subroutines,  5  include  files 
(same  as  UPDATE  common  blocks),  numerous  minor  subroutines.  The  sections  below 
contain  the  descriptions  of  the  following  items: 


•  the  structure  of  the  computer  code, 

•  the  major  subroutines, 

•  the  internal  field  variable  arrays  and  include  file  variables, 

•  the  input/output  files,  and 

•  the  input  variables. 


5.1  Descriptions  of  Subroutines 

Figure  6  shows  the  hierarchy  of  the  subroutines  that  make  up  DUMPSTER.  In  Fig¬ 
ure  7  a  flowchart  shows  the  overall  calling  sequence  of  the  subroutines.  A  brief 
description  of  each  subroutine  is  given  below.  The  intent  of  this  section  is  to  help 
the  user  that  is  interested  in  modifying  DUMPSTER  to  become  acquainted  with  the 
basic  function  of  each  routine.  More  information  about  the  routine  can  be  found  in 
the  FORTRAN77  listing. 

ARRAYS .  Pointer  variables  and  array  sizes  are  calculated  for  the  field  variable 


arrays  (X,  U,  Q,  R,  DU,  DT,  and  TG).  This  subroutine  also  calcu¬ 
lates  the  I  and  J  dimensions  of  each  zone  and  checks  that  the  total 
number  of  computational  cells  does  not  exceed  the  maximum  allowed 
(MAXCELLS). 

BC .  This  subroutine  calls  subroutine  BCZ  for  each  zone. 

BCBWD .  This  subroutine  transfers  the  residual  variable  (R)  from  the  interior 

cells  at  the  bottom  of  a  zone  to  the  boundary  cells  of  the  zone  below 
along  interzone  boundaries  for  the  backward  sweep  of  the  LU-SCJS 


solution. 

BCC .  This  subroutine  calls  subroutine  BCCZ  for  each  zone  to  set  interzone 

boundary  conditions. 
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Figure  6:  Hierarchy  of  DUMPSTER  Subroutines. 
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End 


Figure  7:  A  Flowchart  of  the  Major  Subroutines  of  DUMPSTER. 
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BCCZ 


This  subroutine  transfers  field  variable  data  from  within  one  zone 
to  the  overlaying  boundary  cells  of  a  neighboring  zone  and  also  in 
the  reverse  direction  from  the  interior  cells  of  the  neighboring  zone 
to  the  boundary  cells  of  the  zone.  This  transfer  takes  place  only  on 
boundaries  that  are  specified  as  interzone  boundaries. 

BCFWD .  This  subroutine  transfers  the  residual  variable  (R)  from  the  top  in¬ 

terior  cells  of  a  zone  to  the  bottom  boundary  cells  of  a  neighboring 
zone  above  it  along  interzone  boundaries  for  the  forward  sweep  of 
the  LU-SGS  solution. 

BCZ .  Subroutine  BC  calls  this  routine  for  each  zone  for  setting  boundary 

conditions  in  the  boundary  cells.  The  real  boundary  conditions  that 
affect  the  calculation  of  the  fluxes  are  calculated  in  the  FLXBC  subrou¬ 
tine.  These  boundary  conditions  are  for  plotting  purposes  only  and 
for  transferring  thermodynamic  properties  for  real  gas  calculations. 

CONV .  This  subroutine  calls  subroutine  CONVZ  for  each  zone  to  calculate  the 

level  of  convergence  in  each  zone.  The  zonal  convergence  level  (e„) 
is  defined  as  the  log10  of  the  average  of  the  absolute  value  of  the 
residuals  of  the  mass  density  variable. 

The  average  convergence  level  (CONVA)  is  simply  the  average  of  the 
zonal  convergence  levels  for  all  of  the  zones. 

CONVA  =  i  E 

n=l,N 


N  is  the  total  number  of  zones.  This  subroutine  also  calculates  the 
difference  in  the  total  mass  flow  rate  at  the  left  and  right  boundaries 
of  each  zone  separately.  This  variable  is  only  meaningful  when  a  zone 
is  not  coupled  to  another  zone  and  there  is  inflow  at  the  left  boundary 
and  outflow  at  the  right  boundary  as  for  simple  nozzle  flows. 

CONVZ .  This  subroutine  calculates  the  convergence  level  in  a  zone.  See  des- 

cussion  for  CONV  for  details. 

ELE .  This  subroutine  calculates  coefficients  necessary  to  define  the  stretch¬ 


ing  of  cells  along  a  mesh  line  that  starts  with  exponential  stretching 
then  has  constant  spacing  of  cells  and  then  ends  with  exponential 
stretching.  Typically  this  is  a  stretching  that  might  be  between  two 
flat  plates  having  a  boundary  layer  on  both  plates.  Position  fine  grid 
on  each  boundary  that  stretches  out  to  a  coarser  uniform  grid  in  the 
channel.  The  input  to  this  routine  is:  distance  coordinate  at  the 
bottom  (  Y(I,3)  ),  distance  coordinate  at  the  top  (  Y(I,JDM)  ),  cell 
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spacing  at  the  bottom  (DELYB),  cell  spacing  at  the  top  (DELYT),  and 
the  number  of  cells  in  the  bottom  exponential,  uniform  middle,  and 
top  exponential  regions. 


ERRCHK .  Numerous  checks  are  made  for  the  accuracy  of  the  input  data  in  this 

subroutine.  If  an  error  is  detected  then  a  message  is  printed  to  the 
OPRNT  file  and  execution  is  stopped. 

FLXRI .  This  subroutine  calls  FLXRIZ  subroutine  for  eacli  zone  to  calculate  the 

fluxes  at  the  I-faces  of  all  of  the  cells  and  calculate  the  contribution 
of  the  I-fluxes  to  the  explicit  change  of  the  field  varibles  (DU)  using 
Roe’s  flux-difference  flux  function. 

FLXRIZ .  This  subroutine  calculates  the  fluxes  at  the  I-faces  of  all  of  the  cells 

within  a  zone  using  Roe’s  flux-differencing.  Both  the  inviscid  and 
viscous  fluxes  are  calculated.  The  change  in  the  conservative  solution 
vector  (U)  is  updated  for  the  I-fluxes.  The  fluxes  at  solid  boundaries 
are  directly  set  in  this  routine. 

FLXRJ .  Same  as  FLXRI,  but  for  the  J-faces  of  the  cells. 

FLXRJZ .  Same  as  FLXRIZ,  but  for  the  J-faces  of  the  cells. 

FLXSWI .  Same  as  FLXRI,  except  that  Steger- Warming  flux- vector  splitting  is 

used  to  calculate  fluxes. 

FLXSWIZ .  Same  as  FLXRIZ,  except  that  Steger  Warming  flux-vector  splitting  is 

used  to  calculate  fluxes. 

FLXSWJ .  Same  as  FLXRJ,  except  that  Steger- Warming  flux-vector  splitting  is 

used  to  calculate  fluxes. 

FLXSWJZ .  Same  as  FLXRJZ,  except  that  Steger- Warming  flux-vector  splitting  is 

used  to  calculate  fluxes. 

FLXYI .  Same  as  FLXRI,  except  that  Harten-Yee  flux-vector  splitting  is  used 

to  calculate  fluxes. 

FLXYIZ .  Same  as  FLXRIZ,  except  that  Harten  Yee  flux-vector  splitting  is  used 

to  calculate  fluxes. 

FLXYJ .  Same  as  FLXRJ,  except  that  Harten  Yee  flux-vector  splitting  is  used 

to  calculate  fluxes. 

FLXYJZ .  Same  as  FLXRJZ,  except  that  Harten  Yee  flux- vector  splitting  is  used 

to  calculate  fluxes. 

GSEC .  This  subroutine  is  used  to  call  a  system  function  to  get  the  number 

of  CPU  seconds  that  have  been  used  so  far  in  the  job.  This  subrou¬ 
tine  must  be  customized  for  each  computer  system.  The  standard 
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distribution  subroutine  does  not  call  a  system  routine  and  returns 
0  seconds  of  CPU  time  used.  The  user  inputs  the  maximum  number 
of  CPU  seconds  that  the  job  is  allowed  to  run.  DUMPSTER  period¬ 
ically  checks  the  amount  of  CPU  time  that  remains.  When  the  CPU 
time  is  almost  used  up  DUMPSTER  saves  the  current  solution  and 
terminates  smoothly  so  that  a  restart  run  can  easily  be  made. 


INITU .  This  subroutine  controls  the  initialization  of  the  flowfield.  For  “NEW” 

runs  for  each  zone  this  routine  calls  either  subroutine  UNFICZ,  ONEDZ, 
or  USERICZ  to  initialize  the  flow  field. 

INPUT .  This  subroutine  controls  the  reading  of  the  input  variables  from  the 


input  file  called  IDATA.  An  extended  NAMELIST  type  input  proce¬ 
dure  is  used  to  read  the  data.  The  NAMELIST  routines  are  con¬ 
tained  within  the  two  files  called  XNUTIL  and  XNREAD.  Most  variables 
are  first  set  to  the  default  value.  If  a  variable  is  in  the  input  file, 
then  the  default  value  is  overwritten.  The  value  that  is  read  from 
the  input  file  is  checked  for  valid  range.  Error  messages  are  printed  to 
OPRNT  and  execution  stops  if  input  errors  are  detected  while  reading 
input  data. 


LIMP .  This  subroutine  performs  the  forward  sweep  of  the  LU-SGS  iteration 

at  each  cell  in  a  zone  when  the  k  -  e  model  is  not  activated. 

LIMPKE .  'Phis  subroutine  performs  the  forward  sweep  of  the  LU-SGS  iteration 

at  each  cell  in  a  zone  when  the  k- e  model  is  activated. 

LUDRV .  This  is  the  controller  subroutine  that  calls  the  LIMP,  UIMP,  LIMPKE, 

UIMPKE,  BCFWD,  and  BCBWD  subroutines  that  actually  do  the  LU-SGS 
calculations  within  specific  zones. 

MAIN .  This  is  the  program  module  where  execution  begins.  There  are  sec¬ 


tions  in  this  program  module  that  are  optionally  commented  out  (C 
in  column  1)  that  deal  with  setting  up  the  memory  for  storing  the 
field  variables.  For  most  computer  systems  with  virtual  memory  it 
is  satisfactory  to  hardwire  the  maximum  number  of  cells  (MXNCLS) 
in  a  parameter  statement.  On  Cray  computers  memory  may  be  al¬ 
located  dynamically.  DUMPSTER  will  compute  and  then  request 
the  necessary  storage  for  the  field  variables  automatically.  So  a  user 
only  has  to  compile  the  program  once.  The  dimensions  of  the  field 
variables  are  defined  in  this  routine  either  explicitely  or  dynamically 
depending  upon  which  lines  of  coding  the  user  has  commented  out. 


METR .  This  subroutine  calls  the  subroutine  METRZ  and  METRB  to  calculate 

the  metric  variables  for  each  cell  in  each  zone. 

METRB .  After  the  METRZ  subroutine  calculates  the  metrics  for  each  zone,  this 

subroutine  transfers  metrics  from  the  interior  of  zones  to  the  bound¬ 
ary  cells  of  neighboring  zones  across  interzone  boundaries. 
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METRZ .  This  subroutine  calculates  the  metrics  for  eacli  cell  within  a  zone. 

These  metrics  are  the  face  area  projections  in  the  x  and  r  directions 


for  the  top  and  right  faces,  the  volume  of  the  cell,  and  a  metric  used 
for  calculating  the  axisymmetic  source  term. 

MSHG .  This  subroutine  calls  MSHGZ  to  calculate  a  computational  mesh  in 

each  zone. 

MSHGZ .  This  subroutine  calculates  a  computational  mesh  for  each  zone.  The 

mesh  generation  procedure  is  described  in  Section  5.10. 

ONEDZ .  This  subroutine  calculates  initial  conditions  for  a  zone  using  one¬ 


dimensional  gas  dynamics  relations.  These  initial  conditions  are  use¬ 
ful  for  setting  initial  conditions  for  nozzles  and  combustors  with  sim¬ 
ple  geometry. 


PBNR .  This  utility  routine  prints  a  openning  banner  to  the  print  out  fde 

OPRNT. 

PHDR .  This  utility  routine  advances  the  print  out  to  a  new  page  and  writes 

the  TITLE  character  string  at  the  top  of  the  page. 

PLOT .  This  subroutine  writes  the  header  information  to  the  plot  file  (OPLTF) 


in  the  format  for  TECPLOT  (Amtec  Engineering’s  interactive  plot¬ 
ting  program.)  PLOT  then  calls  PLOTZ  for  each  zone  to  write  out  the 
field  variable  data  to  the  plot  file.  PLOT  is  always  called  at  the  end  of 
a  calculation.  The  user  may  specify  in  the  input  data  to  create  plot 
files  at  intermediate  steps  during  the  calculation.  These  plots  are 
used  for  plotting  meshes,  contours,  and  vector  fields  for  visualizing 
the  calculated  solution  data. 

PLOTH .  During  the  calculation  DUMPSTER  may  (depending  upon  what  the 

user  specifies  in  the  input  data)  write  out  the  TAP  data  to  a  scratch 
file  from  subroutine  WTAP.  At  the  end  of  the  calculation  the  PLOTH 
subroutine  is  called  to  read  the  TAP  data  from  the  scratch  file,  reor¬ 
ganize  it,  and  then  write  it  out  to  the  file  called  OPLTH  in  TECPLOT 
format  ready  for  plotting.  TAP  data  consist  of  the  field  variables 
at  a  set  of  mesh  cells  at  various  computational  steps.  These  data  is 
useful  for  plotting  convergence  history  of  various  field  variables. 


PLOTZ .  This  subroutine  writes  out  the  field  variables  data  for  a  zone  for 

creating  field  plots. 

PMJR .  This  subroutine  controls  the  printing  of  the  major  printout.  It  calls 

PMJRZ  for  each  zone  to  print  out  field  variable  data. 

PMJRZ .  At  user  specified  columns  (I=constant)  for  a  zone  the  coordinates,  ve¬ 

locity  components,  pressure,  temperature,  density,  turbulence  quan¬ 
tities,  total  pressure,  and  total  temperature  are  printed  to  file  OPRNT 
for  each  cell. 
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PMNR 


This  subroutine  controls  the  printing  of  the  minor  print  out.  This 
is  a  selected  amount  of  data  that  is  printed  several  time  during  a 
calculation  to  monitor  the  solution  progress.  PMNRZ  is  called  for  each 
zone. 


PMNRZ .  This  subroutine  prints  various  field  variables  for  a  zone  along  th< 

I-direction  of  the  zone  on  the  top,  middle,  and  bottom  boundaries. 
Field  variables  are  also  printed  at  user  selected  I-columns. 

PPRP .  This  subroutine  controls  the  printing  of  thermodynamic  and  trans¬ 

port  properties.  PPRPZ  is  called  for  each  zone. 

PPRPZ .  The  thermodynamic  and  transport  properties  are  printed  to  OPRNT 

for  a  zone  at  the  I-columns  specifed  for  minor  prints  in  the  input 
data. 

PQUIK .  This  subroutine  controls  the  printing  of  a  quick  and  short  amount 


of  data  to  the  OPRNT  file  and  to  standard  out  for  monitoring  the 
progress  of  the  solution  procedure.  The  computational  step,  CFL 
mulitplier,  and  the  convergence  level  is  printed.  PQUIKZ  is  called  for 
each  zone  to  print  out  a  sampling  of  the  flow  field  data. 


PQUIKZ .  This  subroutine  writes  out  a  sampling  of  the  data  for  a  zone.  This 

subroutine  is  easily  modifiable  so  that  the  user  may  customize  it  for 
various  calculations. 

PREP .  After  the  input  data  has  been  read,  and  the  array  memory  allocated, 

the  program  calles  PREP  subroutine  to  calculate  and  initialize  various 
parameters  and  variables.  These  include  setting  metrics  between 


zones  along  interzone  boundaries,  setting  up  temperature- dependent 
property  coefficients,  writing  header  information  to  plot  files,  and 
initializing  the  step  counter  and  the  CFL  multiplier. 

PROC .  The  main  iteration  loop  in  which  ISTEP  in  the  index  variable  is  con¬ 

tained  in  this  subroutine.  Each  pass  through  the  loop  constitutes  an 
iteration  or  computational  step.  The  CFL  multiplier  and  the  relax¬ 
ation  coefficient  are  automatically  adjusted  in  this  routine  to  help 
maintain  stability.  If  the  time  step  falls  below  the  explicit  stability 
limit  the  code  will  deactivate  the  implicit  solution  procedure 


PXY .  The  user  has  the  option  to  print  out  the  x  and  y  (r)  coordinates  at 

the  beginning  of  the  calculation.  This  subroutine  calls  the  subroutine 
PXYZ  to  print  out  the  x  and  y  coordinates  for  each  zone. 

PXYZ .  The  x  and  y  coordinates  for  a  zone  are  printed  out  to  file  OPRNT  in 

this  subroutine. 

RMSH .  The  computational  mesh  may  be  read  from  an  input  file  IMESH  or 

generated  internally.  RMSH  is  the  subroutine  that  controls  the  reading 
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of  the  external  mesh.  It  calls  subroutine  RMSHZ  to  read  the  mesh  for 
each  zone. 


RMSHZ .  This  subroutine  reads  the  mesh  coordinates  for  a  zone  from  the  input 

file  IMESH. 

RSTR .  At  the  end  of  a  calculation  DUMPSTER  saves  the  state  of  the  flow 


field  in  a  restart  file  called  ORSTR.  This  subroutine  can  read  in  this 
restart  file  (renamed  as  IRSTR)  to  continue  with  the  calculation  at 
a  later  time.  The  restart  file  contains  the  field  variables  X  and  U 
and  some  information  about  the  calculation  when  it  was  saved  that 
will  allow  DUMPSTER  to  continue  the  calculation  as  if  it  were  not 
interupted. 

SETBCL .  DUMPSTER  calls  SETBCL  at  the  beginning  of  the  calculation  to  set 

up  some  arrays  that  store  the  boundary  conditions  at  the  left  inflow 
boundary  (1=1).  For  subsonic  inflow  boundary  conditions  this  in¬ 
cludes  flow  angle,  total  pressure,  and  total  temperature  distributions 
in  the  radial  direction.  In  this  subroutine  these  input  arrays  are  in¬ 
terpolated  into  another  set  of  arrays  that  have  a  one  to  one  mapping 
to  the  boundary  cells  at  the  inflow  boundary.  When  boundary  con¬ 
ditions  are  applied  at  the  inflow  boundary,  the  flow  angle  and  total 
conditions  are  readily  available  for  each  cell. 


SETBCLZ .  This  subroutine  calculates  the  boundary  conditions  arrays  for  the 

left  inflow  boundaries  for  a  zone. 

SOLVE .  This  subroutine  is  a  utility  subroutine  that  solves  a  system  of  N  linear 

algebraic  equations  using  a  form  of  Gaussian  elimination. 

STOR .  At  the  end  of  a  calculation  this  subroutine  stores  the  flow  field  in  a 

restart  file  called  ORSTR. 

TBLU .  This  utility  function  performs  linear  and/or  quadratic  interpolation 

into  a  table  of  floating-point  numbers. 

TBLUI .  This  utility  function  is  the  same  as  TBLU  except  that  it  interpolates 

into  a  table  of  two  variables — one  integer  and  one  floating-  point. 

TKE .  This  subroutine  controls  the  calls  to  subroutine  TKEZ  to  calculate  the 

eddy  viscosities  when  using  the  fc-c  turbulence  model. 

TKEZ .  This  subroutine  calculates  the  eddy  viscosity  from  the  values  of  the 

turbulent  kinetic  energy  (k)  and  turbulence  dissipation  rate  (c)  in  a 
zone  when  the  fc-e  turbulence  model  is  used.  The  values  of  k  and  t 
are  also  calculated  near  walls  for  coupling  with  the  wall  function. 

TSSET .  When  the  time  step  is  uniform,  this  subroutine  is  used  to  set  the 

time  step  in  a  zone.  It  is  called  from  subroutine  TSTP. 
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TSTP 


This  subroutine  is  used  to  control  the  calculation  of  the  time  step  at 
each  computational  cell.  The  time  step  may  be  uniform  or  based  on 
local  conditions. 

TSTPZ .  This  subroutine  calculates  the  local  time  step  in  each  cell  of  a  zone. 

TURB .  This  subroutine  controls  the  calculation  of  the  eddy  viscosity  in  the 

zones  when  using  the  Baldwin-Lomax  turbulence  model,  it  calls 
TURBZ  for  each  zone. 

TURBZ .  This  subroutine  calculates  the  eddy  viscosity  and  eddy  thermal  con¬ 

ductivity  in  a  zone.  Only  boundary  layers  on  upper  (J=JD)  and 
lower  (J  =  l)  walls  are  allowed. 

UIMP .  This  subroutine  performs  the  backward  sweep  of  the  LU-SGS  itera¬ 

tion  at  each  cell  in  a  zone  when  the  k-e  model  is  not  activated. 

UNFICZ .  This  subroutine  sets  the  field  variables  in  a  zone  to  some  uniform 

initial  conditions. 

UPDQ .  This  subroutine  controls  the  calls  to  subroutine  UPDQZ  to  sets  the 

thermodynamic  and  transport  properties. 

UPDQZ .  This  subroutine  calculates  and  set  the  thermodynamic  and  transport 

properties  in  a  zone. 

UP  DU .  This  subroutine  controls  the  update  of  the  conservative  field  variable 

array  U  at  the  end  of  a  computational  step.  Subroutine  UPDUZ  is 
called  to  update  the  U-array  for  each  zone. 

UPDUZ .  The  conservative  field  variable  array  U  is  updated  in  this  subroutine. 

DU  is  added  to  the  U-array  from  the  last  time  step.  This  subroutine 
also  may  limit  the  size  of  the  change  in  U  (i.e.  DU)  to  a  percentage 
of  U. 

USERICZ .  The  user  may  optionally  add  custom  coding  to  DUMPSTER  for  ini¬ 

tializing  the  flow  field.  This  coding  may  be  added  to  this  subroutine. 

WERRQR .  This  utility  subroutine  writes  error  messages  to  the  print  file  OPRNT. 

WMSSG .  This  utility  subroutine  writes  messages  to  the  print  file  OPRNT. 

WTAP .  This  subroutine  controls  the  writting  of  tap  data  to  the  output  file 

OTAP.  It  calls  the  subroutine  WTAPZ  to  write  out  the  data  for  each  tap 
in  each  zone. 

WTAPZ .  When  this  subroutine  is  called  it  writes  out  the  field  variables  at 

each  tap  for  a  zone.  A  tap  is  just  a  mesh  cell  location  where  the  user 
may  request  (via  the  input  data  file  IDATA)  that  field  variable  data 
be  stored  periodically  during  a  calculation.  After  the  calculation  is 
completed,  this  data  is  organized  for  plotting  with  TECPLOT  and 
written  to  the  plot  file  OPLTH. 
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WWARNG .  This  subroutine  writes  out  a  warning  message  to  the  print  file  QPRNT. 

XINT2 .  This  subroutine  is  used  to  interpolated  from  a  non-rectangular  2D 

mesh  to  a  point.  It  returns  the  cell  where  the  interpolation  was  done 
and  the  weighting  coefficients  for  the  bi-linear  interpolation.  These 
coefficients  are  used  in  the  calling  subroutine  XPATCH  to  calculate  the 
values  of  the  field  variables  at  the  interpolated  point. 

XINT2S .  This  subroutine  sets  up  the  arrays  for  XINT2. 

XPATCH .  DUMPSTER  may  restart  a  calculation  by  initializing  the  flow  field 

by  interpolation  from  a  flow  field  on  a  coarser  mesh.  This  process  is 
controlled  in  this  subroutine.  It  calls  XINT2S  and  XINT2  to  perform 
the  searching  and  interpolations. 

YGEOM .  The  user  may  supply  custom  coding  for  generating  complex  geome¬ 

tries.  This  coding  should  be  placed  in  this  subroutine. 

YMSH .  This  subroutine  calculates  a  distribution  of  mesh  points  along  a  ver¬ 

tical  line.  It  is  called  from  MSHGZ. 

(XN  ROUTINES) .  The  extended  NAMELIST  routines  are  contained  in  the  two  files: 

XNREAD  and  XNUTIL.  These  routines  should  not  be  altered. 

ZERO .  This  routine  initializes  the  arrays  to  zero  at  the  start  of  a  calculation. 
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5.2  Description  of  Field  Variable  Arrays 


This  section  gives  a  brief  description  of  the  field  varibles  that  are  used  in  the  DUMP¬ 
STER  coding.  Field  variables  are  those  that  are  defined  at  each  computational  cell. 


Internal 
Array  Name 

Math 

Symbol 

Description 

X 

X 

Metrics:  coordinates,  area  projections,  and 

volumes 

U 

u 

Conservative  fluid  dynamic  variables 

DU 

A  U 

Changes  in  conservative  fluid  dynamic  variables. 

R 

R 

Residual  of  fluid  dynamics  equations  (similar  to 
DU) 

Q 

Q 

Thermodynamic  and  transport  properties 

DT 

St 

Local  time  steps.  (One  element.) 

TG 

Gk 

Turbulent  dissipation  function  used  in  the  k-e 
model.  (One  dimensional) 

Each  of  the  above  arrays  contain  many  elements  as  shown  in  the  following  list.  The 
last  subscript  in  the  arrays  shown  below  are  integer  parameters  defined  in  Section  5.4. 


X(I ,  J  ,NX) .  (a;)  Axial  coordinate 

X(I,J,NY) .  (r  or  y)  Radial  or  vertical  coordinate 

X(I,J.NAXT) .  (Aj)  Area  projection  of  the  face  at  the  top  i,j  +  ^  of  cell  i,j  in  the 

x  direction 

X(I,J,NAYT) .  (Ay  )  Area  projection  of  the  face  at  the  top  ( i,j  +  of  cell  i,j  in 

the  y  direction 

X(I,  J.NAXR) .  (A^)  Area  projection  of  the  face  at  the  right  side  (i  +  \,j)  of  cell  i,j 

in  the  x  direction 

X(I, J.NAYR) .  (A^)  Area  projection  of  the  face  at  the  right  side  (t  +  j,j)  of  cell  i,j 

in  the  y  direction 

X(I,J,NV0L) .  (V)  Volume  of  cell  i,j 

X(I,J,NAZ) .  (A*)  Area  projection  of  front  and  back  of  cell  i,j  in  the  y  direction 

U(I ,  J , NR) .  (p)  Mass  density 

U(I,J,NRU) .  (pu)  x  momentum 
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U(I,J,NRV) .  (pv)  y  momentum 

U(I,J,NRW) .  (pu>)  Swirl  momentum 

U(I,  J,NE) .  (£)  Total  energy  density,  E  =  pe  +  |(u2  +  t>2  +  w2) 

U(I,J,MRK) .  ( pk )  Turbulence  energy 

U(I,J,NRE) .  (pe)  Turbulence  dissipation  rate 

U(I,J,NRF1) .  ( pfi )  Species  number  1 

U(I,J,NRF2) .  (ph)  Species  number  2 

U(I,J,NRF3) .  (p/3)  Species  number  3 

U(I,J,NRF4) .  (pf 4)  Species  number  4 

U(I,J,NRF5) .  (p/5)  Species  number  5 

Q(I,  J.NVIS) .  (p)  Effective  viscosity 

Q(I, J.NCND) .  ( k )  Effective  conductivity 

Q(I,J,NGAM) .  (7)  Effective  ratio  of  specific  heat  capacities  for  thermal  equation  of 

state.  7  =  P/(pe) 

Q (I ,  J.NEXP) .  (7)  Defined  as  the  ratio  of  the  speed  of  sound  squared  times  density 

to  the  pressure.  7  =  ^ 

Q(I,J,NCVM) .  (Cv)  Mean  specific  heat  capacity  at  constant  volume 

Q(I,J,NCHI) .  (x)  Thermodynamic  function  used  in  Roe’s  fluxes 

Q(I,J,NKPA) .  (k)  Thermodynamic  function  used  in  Roe’s  fluxes 


Since  the  elements  in  the  arrays  U,  DU,  and  R  correspond  one  to  one,  only  the  U 
elements  of  the  array  are  listed  above.  For  an  example:  U(I,J,NRW)  is  the  swirl 
momentum,  DU(I,J,NRW)  is  the  change  in  the  swirl  momentum,  and  R(I,J,NRW) 
is  the  residual  of  the  swirl  momentum  equation  (almost  the  same  as  DU).  The  third 
array  index  is  specified  indirectly  using  a  variable  such  as  NRW  in  the  above  example. 
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5.3  Description  of  the  AGLBL  Include  File 


The  include  file  AGLBL  is  included  in  all  subroutines  and  functions.  It  contains  the 
following  coding. 


IMPLICIT  REAL*8  (A-H.O-Z) 
COMMON  /DEBUG/  IDBUG(IO) 


This  file  makes  all  floating  point  variables  double  precision  on  a  32-bit  computer.  To 
compile  DUMPSTER  on  a  64-bit  computer  like  the  Cray-Y/MP,  the  above  include 
file  must  be  changed  to  the  following. 


IMPLICIT  REAL  (A-H.O-Z) 
COMMON  /DEBUG/  IDBUG(IO) 


5.4  Description  of  the  AGNRL  Include  File 

AGNRL  contains  numerous  parameters  that  either  do  not  change  during  a  calculation 
or  change  only  rarely.  The  important  parameters  in  this  include  file  are  listed  below. 
The  standard  settings  for  version  1.0  of  DUMPSTER  for  these  variables  are  indicated. 
These  settings  may  be  changed. 


MAXZON .  Maximum  number  of  zones  (MAXZ0N=3) 

MAXZONP .  Maximum  number  of  zones  plus  one  (MAXZONPs^ 

MAXID .  Maximum  number  of  cells  in  i-direction  (MAXID=200) 

MAXJD .  Maximum  number  of  cells  in  j-direction  (MAXJD=50) 

MAXNBC .  Maximum  number  of  points  used  to  set  inflow  boundary  conditions 

distribution  (MAXNBC=100) 

MAXI JD .  Maximum  of  MAXID  and  MAXJD  (MAXIJD=200) 

MAXGEO .  Maximum  of  number  of  geometry  data  points  (MAXGE0*5O) 

MAXCOL .  Maximum  number  of  columns  in  minor  printouts  (MAXCOL“10) 
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MAXPET .  Maximum  number  of  printouts  (MAXPRT=100) 

HAXPLT .  Maximum  number  of  plots  (MAXPLT=100) 

HAXTPS .  Maximum  number  of  field  variable  taps  (MAXTPS*20) 

MAXNEX .  Maximum  number  of  elements  in  X-array  (MAXNEX=8) 

MAXNEU .  Maximum  number  of  elements  in  U-array  (MAXNEU*8) 

MAXNEQ .  Maximum  number  of  elements  in  Q-array  (MAXNEQ=8) 

MAXNEDT .  Maximum  number  of  elements  in  DT-array  (MAXNEDT=1) 

MAXNETG .  Maximum  number  of  elements  in  TG-array  (MAXNETG=1) 

MAXNSPC .  Maximum  number  of  species  in  U-array  (MAXNSP05) 

MAXNKEU .  Maximum  number  of  elements  used  for  the  k  t  turbulence  model 

(MAXNKEU=2) 

MAXCLLS .  Maximum  number  of  cells  as  set  in  subroutine  MAIN 

IDMX .  Dimension  for  X-array  as  set  in  subroutine  ARRAYS 

IDMU .  Dimension  for  U-array  as  set  in  subroutine  ARRAYS 

IDMq .  Dimension  for  Q-array  as  set  in  subroutine  ARRAYS 

IDMR .  Dimension  for  R- array  as  set  in  subroutine  ARRAYS 

IDMDU .  Dimension  for  DU-array  as  set  in  subroutine  ARRAYS 

IDMDT .  Dimension  for  DT-array  as  set  in  subroutine  ARRAYS 

IDMTG .  Dimension  for  TG-array  as  set  in  subroutine  ARRAYS 

AXI .  Coordinate  switch:  =1  for  Axisymmetric,  =0  for  2D  planar 

TWD .  Coordinate  switch:  =0  for  Axisymmetric,  =1  for  2D  planar 

LOSTD .  Standard  out  file  (L0STD=6) 

LOSCRA .  Scratch  file  (L0SCRA=9) 

LIDATA .  IDATA  input  file  (LIDATA=10) 

LIMESH .  IMESH  input  file  (LIMESH=11) 

LIRSTR .  IRSTR  input  file  (LIRSTRM2) 

LORSTR .  ORSTR  output  file  (L0RSTR=2O) 

LOPRNT .  OPRNT  output  file  (L0PRNT=21) 
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LOPLTF .  OPLTF  output  file  (LOPLTF-22) 

LOPLTH .  OPLTH  output  file  (L0PLTH*23) 

LOTAPS .  OTAPS  output  file  (LOTAPS-30) 

LOCONV .  OCONV  output  file  (L0C0NV-31) 


The  following  parameters  in  the  AGNRL  include  file  axe  used  as  subscripts  in  the 
field  variable  arrays.  For  example  the  x-component  of  momemtum  in  the  U-array  is 
U(I,J,NRU)  whose  mathematical  symbol  is  pu. 


NX .  x-coordinate  (NX  =  1) 

NY .  y-  or  r-coordinate  (NY=2) 

NAXT .  Area  Projection  in  x-direction  at  top  face  (NAXT*3) 

NAYT .  Area  Projection  in  y-direction  at  top  face  (NAYT=4) 

NAXR .  Area  projection  in  x-direction  at  right  face  (NAXR*5) 

NAYR .  Area  projection  in  y-direction  at  right  face  (NAYR*6) 

NVOL .  Volume  of  cell  (NV0L=7) 

NAZ .  Area  projection  of  sides  of  cell  in  y-direction  (NAZ*8) 

NR .  Mass  density  (NR*  1) 

NRU .  x-momentum  (NRU=2) 

NRV .  y-momentum  (NRV=3) 

NRW .  Swirl  momentum  (NRW=4) 

NE .  Total  energy  densiy  (NE=5) 

NRK .  Turbulence  kinetic  energy  pk  (NRK=6) 

NRE .  Turbulence  dissipation  rate  pi  (NRE*7) 

NRF1 .  Mass  fraction  for  species  1  (NRF1»8) 

NRF2 .  Mass  fraction  for  species  2  (NRF2*9) 

NRF3 .  Mass  fraction  for  species  3  (NRF3*10) 

NRF4 .  Mass  fraction  for  species  4  (NRF4=11) 
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HRF5 .  Mass  fraction  for  species  5  (NRF5*12) 

NVIS .  Viscosity  (NVIS*1) 

HCMD .  Thermal  conductivity  (NCND*2) 

NGAM .  Ratio  of  specific  heat  capacities  =  1  +  £  (NGAM*3) 

NCVN .  Mean  specfic  heat  capacity  at  constant  volume  (NCVM*4) 

NEXP .  Ratio  of  pc J  to  P  (NEXP-5) 

NfCHI .  Thermodynamic  function  used  for  Roe’s  fluxes  (HCHI*6) 

BKPA .  Thermodynamic  function  used  for  Roe’s  fluxes  (NKPA*7) 


5.5  Description  of  the  GCTRL  Include  File  Variables 

The  variables  in  the  GCTRL  include  file  are  described  below.  These  variables  are  global 
in  use  ( i.e .  they  are  applicable  to  all  zones.) 


ICHTRL .  Execution  Method 

NSTEPS .  Number  of  steps  to  take  for  this  run 

NZONES .  Total  number  of  zones  that  are  used 

ICRD .  Coordinate  system  switch:  =0  for  2D  planar,  =1  for  axisymmetric 

NEU .  Number  of  elements  in  the  U-array.  Normally  NEU=5  for  p,pu,  pv,  pw, 

and  E 

C NVTL .  Convergence  tolerance.  The  number  of  orders  of  magnitude  that  the 

convergence  level  must  be  reduced  before  convergence  is  reached 

SECJOB .  Alloted  number  of  seconds  for  the  job  to  run 

ITHRRO .  Thermodynamic  model  selection  integer 

IVISC .  Viscosity  model  selection  integer 

ITHINL .  Thin-layer  option  for  diffustion  terms 

ICOND .  Thermal  conductivity  model  selection  integer 

NSPC .  Number  of  chemical  species  being  transported 

NTKE .  Number  of  field  variables  used  in  turbulence  model.  This  is  normally 

0  or  2  (for  k-(  model) 


GAMC .  Ratio  of  specific  heat  capacities 

GASC .  Perfect  gas  constant  (Normally  =  287.Joules/Kg-Deg.K) 

PRDL .  Laminar  Prandtl  number 

PRDT .  Turbulent  Prandtl  number 

CVC0N(4) .  Coefficients  for  the  specific  heat  capacity 

CVMC0NC4) .  Coefficients  for  the  mean  specific  heat  capacity 

VISC0N(6) .  Coefficients  for  the  viscosity 

CNDC0N(6) .  Coefficients  for  the  thermal  conductivity 

SPCCON(MAXNSPC)  . . .  Diffusion  coefficients  for  species  transport 

TKEC0N(5) .  Coefficients  for  the  k-(  turbulence  model 

ICDT .  Control  integer  for  time  step  calculation 

IFLXF .  Flux  function  control  interger 

RFGS .  Gauss- Seidel  relaxation  coefficient 

RFTS .  Time  step  relaxation  coefficient 

BETA .  Coefficient  in  the  LU-SGS  relaxation 

CHGMAX .  Fraction  of  change  allowed  when  updating  U-array 

CFLBEG .  CFLM  used  at  the  beginning  of  the  calculation 

CFLFAC .  A  factor  that  multiplies  the  CFLM  to  increase  it  from  one  step  to 

the  next 

CFLMAX .  The  maximum  value  for  the  CFLM 

FLXCFS .  Coefficients  used  in  the  flux  functions 

IMSHG .  Mesh  generation  method  control  integer 

ISTEP .  Current  computational  step  counter 

ICTOL .  Control  integer  for 

IMPL .  Impicit/Explicit  switch:  =0  for  explicit  =1  for  implicit 

CONVA .  Current  average  convergence  level 

CFLM .  Current  CFL  multiplier 

TIME .  Current  time 


53 


SECCPU 


Current  number  of  seconds  used  in  calculation 


NPLT .  N  umber  of  plots  to  be  make 

IPRT .  Number  of  minor  print  outs  to  be  make 

NPRT .  Steps  at  which  minor  print  outs  are  to  be  made 

IPLT .  Steps  at  which  plots  are  to  be  made 

IFQPX .  Interval  of  steps  for  quick  print 

ISKP .  Skip  interval  for  the  i-index  for  major  print  outs 

JSKP .  Skip  interval  for  the  ./-index  for  major  print  outs 

IPRTQ .  Properties  print  out  control  integer 

IPRTM .  Mesh  coordinates  and  metrics  print  out  control  integer 

5.6  Description  of  the  ZCTRL  Include  File  Variables 

ZCTRL  include  file  contains  variables  that  define  attributes  of  each  zone.  The  di¬ 
mension  of  the  variables  are  also  shown.  MAXZON  is  the  maximum  number  of  zones 
allowed  and  MAXZONP  is  equal  to  MAXZON  plus  one. 

IDZ(MAXZON) .  Number  of  cells  in  i-direction  (including  four  boundary  cells) 

JDZ(MAXZON) .  Number  of  cells  in  ./-direction  (including  four  boudary  cells) 

IPNTX (MAXZONP) _  Array  pointer  variable  for  X  array 

IPNTU(MAXZONP) _  Array  pointer  variable  for  U-array 

IPNTQ (MAXZONP) _  Array  pointer  variable  for  Q-array 

IPNTR (MAXZONP) _  Array  pointer  variable  for  R-array 

IPNTDU (MAXZONP) . . .  Array  pointer  variable  for  DU-array 
IPNTDT (MAXZONP) . . .  Array  pointer  variable  for  DT-array 
IPNTTG(MAXZONP)  . . .  Array  pointer  variable  for  TG-array 

ISHUMZ(3  .MAXZON)  .  Difference  in  i  index  for  neighboring  cells  that  are  on  opposite  sides  of 
an  interzone  boundary  Upper  zone  i-index  minus  lower  zone  i-index. 

IIC(MAXZON) .  Initial  conditions  control  integer 


PRESI (MAXZON) .  Initial  condition  pressure 

TEMPI (MAXZON) .  Initial  condition  temperature 

VVELI(3 .MAXZON) .. .  Initial  condition  velcity  components 

TKI  (MAXZON) .  Initial  condition  turbulence  kinetic  energy 

TEI  (MAXZON) .  Initial  condition  turbulence  dissipation  rate  energy 

FMI(MAXNSPC, MAXZON)  Initial  condition  species  mass  fraction 

THMACH (MAXZON) ....  Throat  Mach  number  sometimes  used  for  initial  conditions 

KBCB( 3, MAXZON) _  Boundary  condition  type  on  bottom  boundary  segments 

KBCT (3, MAXZON)  ....  Boundary  condition  type  on  top  boundary  segments 


KBCL (MAXZON) .  Boundary  condition  type  on  left  boundary 

KBCR(MAXZON) .  Boundary  condition  type  on  right  boundary 


IBCB (4, MAXZON)  ... .  Cell  index  array  for  defining  start  and  end  of  boundary  condition 
segments  on  the  bottom  boundary 

IBCT (4 .MAXZON) -  Cell  index  array  for  defining  start  and  end  of  boundary  condition 

segments  on  the  top  boundary 

NTRPBC (MAXZON)  ....  Interpolation  option  for  setting  up  the  boundary  conditions  at  the 
left  inflow  boundary 

NBCVPT(MAXZON)  ....  The  number  of  points  that  define  the  boundary  condition  array  for 
the  left  boundary  for  the  velocity  angles,  pressure,  and  temperature 

NBCTKE (MAXZON) ....  The  number  of  points  that  define  the  boundary  condition  array  for 

the  left  boudnary  for  the  turbulence  model 

NBCSPC (MAXZON) ....  The  number  of  points  that  define  the  boundary  condition  array  for 
the  left  boundary  for  the  species  mass  fraction 

PEX(MAXZON) .  The  static  pressure  for  outflow  boundary  conditions 

TEX(MAXZON) .  Static  temperature  for  outflow  boundary  conditions 

BCCFS( 9, MAXZON)  . . .  Boundary  condition  coefficients 

TWL(MAXZON) .  Temperature  of  left  wall 

TWR(MAXZON) .  Temperature  of  right  wall 

TWB(MAXZON) .  Temperature  of  bottom  wall 

TWT(MAXZON) .  Temperature  of  top  wall 
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BCVPT(MAXNEU+1 .MAXNBC .MAXZON)  Boundary  conditions  array  for  defining  velocity  compo¬ 
nents,  pressure,  and  temperature  distributions  at  the  left  boundary 

BCTKE(MAXNKEU+1, MAXNBC, MAXZON)  Boundary  conditions  array  for  defining  turbulence 
quantities  at  the  left  boundary 

BCSPC(MAXNSPC+1, MAXNBC, MAXZON)  Boundary  conditions  array  for  defining  species  mass 
fractions  at  the  left  boundary 


IDIM(MAXZON) .  Number  of  internal  cells  in  t-direction 

JDIM(MAXZON) .  Number  of  internal  cells  in  j-direction 

NSB(MAXZON) .  Number  of  points  defining  the  geometry  of  the  bottom  surface  of 

zone 


NTRPGB (MAXZON) ....  Interpolation  method  for  bottom  geometry 

NST(MAXZON) .  Number  of  points  defining  the  geometry  of  the  top  surface  of  zone 

NTRPGT (MAXZON)  ....  Interpolation  method  for  top  geometry 

ISN(MAXZON) .  Control  switch  for  generating  mesh  in  the  t-direction 

NCX(4, MAXZON) .  Number  of  cells  in  four  segments  of  mesh  in  the  t-direction 

JSW(MAXZON) .  Control  switch  for  generating  mesh  in  the  7-direction 

NJBK (MAXZON) .  Number  of  blocks  of  cell  in  the  7-direction 

JBK(5, MAXZON) .  Number  of  cells  in  block 

NIJBKMAXZON) .  Number  of  t-indices  for  defining  mesh 

IJBK(10,MAXZ0N)  ...  An  array  of  t-indices  where  mesh  is  defined 

XYCONV(MAXZON) -  Conversion  factor  that  multiplies  coordinates  of  externally  generated 

mesh  when  it  is  read  into  DUMPSTER 

XSBO  (MAXZON) .  x-coordinate  of  the  origin  of  the  geometry  array  for  bottom  surface 

YSBO (MAXZON) .  y-coordinate  of  the  origin  of  the  geometry  array  for  bottom  surface 

XSTO(MAXZON) .  x-coordinate  of  the  origin  of  the  geometry  array  for  top  surface 

YSTO (MAXZON)  ....  . .  y-coordinate  of  the  origin  of  the  geometry  array  for  top  surface 

XSB(MAXGEO  .MAXZON)  x-coordinate  of  the  geometry  array  for  bottom  surface 
YSB(MAXGEO , MAXZON)  y-coordinate  of  the  geometry  array  for  bottom  surface 
XST(MAXGEO  .MAXZON)  x-coordinate  of  the  geometry  array  for  top  surface 


YST(MAXGEO,MAXZON)  y-coordinate  of  the  geometry  array  for  top  surface 


XCTR(HAXZON) .  x-coordinate  of  the  mesh  origin 

DELX(MAXZON) .  Cell  spacing  in  the  x-direction  at  the  mesh  origin 

SLX(4  ,MAXZ0N) .  A  lenth  of  the  mesh  segement  of  a  stretch  factor  for  the  mesh  segment 

in  the  x-direction 


DYBJBK(10, MAXZON)  Cell  spacing  in  the  y-direction  at  the  bottom  of  the  zone 
DYTJBK( 10, MAXZON)  Cell  spacing  in  the  y-direction  at  the  top  of  the  zone 
EJBK(lO,5,MAXZON)  Stretch  factors  for  mesh  generation  in  the  y-direction 


ITCHK .  Check  sum  for  number  of  taps 

NTAPS(MAXZON) .  Number  of  taps  in  zone 


ITAP(MAXTPS .MAXZON)  i-indices  of  taps 

JTAP (MAXTPS , MAXZON )  ji-indices  of  taps 

NCOLPRT (MAXZON)  . . .  Number  of  columns  in  minor  print 

ICOLPRT(MAXCOL, MAXZON)  t-indices  of  columns  in  minor  print 

XMFR(MAXID .MAXZON)  Mass  flow  rate  across  an  t-plane 

CNVP (4, MAXZON)....  Convergence  parameters 


5.7  Input  Files 


DUMPSTER  reads  data  from  one  or  more  of  the  files  listed  in  the  table  below.  When 
running  DUMPSTER  the  files  IRSTR  and  IMESH  must  be  present  even  if  they  are  not 
used. 


Logical  Unit  Name 

Description 

Type 

File  Name 

LIDATA 

Input  Variable  File 

ASCII 

I  DATA 

LIRSTR 

Input  Restart  Data  File 

binary 

IRSTR 

LIMESH 

Input  External  Mesh  File 

ASCII 

IMESH 

5.8  Output  Files 


DUMPSTER  writes  data  to  one  or  more  of  the  files  listed  in  the  table  below. 


Logical  Unit  Name 

Description 

Type 

File  Name 

LOPRNT 

Output  Printer  File 

ASCII 

LORSTR 

Output  Restart  File 

binary 

LOPLTF 

Output  Plot  File 

ASCII 

LOPLTH 

Output  Plot  History  File 

ASCII 

OPLTH 

LOTAPS 

Output  Taps  History  File 

ASCII 

OTAPS 

LOCONV 

Output  Convergence  File 

ASCII 

OCONV 

LOSTD 

Standard  Output 

ASCII 

(Unit  6) 

The  file  OPRNT  is  suitable  for  printing  on  a  line  printer  or  viewing  on  a  computer  screen. 
QRSTR  may  be  renamed  IRSTR  for  restarting  a  calculation.  Files  OPLTH,  OPLTF,  OCONV, 
and  OTAPS  are  in  an  ASCII  format  that  is  readable  by  the  plotting  program  called 
TECPLOT[40]. 

The  logical  unit  name  LOSCRA  is  used  for  a  scratch  file  called  SCRA.  This  file  is  used 
for  temporary  storage  of  the  tap  data. 
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5.9  Dimensional  Units 


The  dimensional  units  of  the  input  and  output  variables,  as  well  as  the  internal  vari¬ 
ables  in  DUMPSTER,  are  in  MKS  (Meter-Kilogram-Second }.  The  units  for  common 
variables  are  tabulated  below. 


Physical  Quantity 

Units  in  MKS 

Basic  Units 

Length 

meters 

m 

Area 

m2 

Volume 

m3 

Time 

seconds 

s 

Velocity 

m/s 

Mass 

kilogram 

Kg 

Density 

I<g/m 3 

Temperature 

deg.  Kelvin 

°K 

Force 

Newtons  (N) 

Kg  •  m/s 2 

Pressure 

Pascals  (Pa) 

N/m 2 

Energy 

Joules  (J) 

N  ■  m 

Energy  density 

J/Kg  —  m2/s2 

Heat  transfer  rate 

J/(m2  •  s) 

Viscosity  (absolute) 

I<g/(m  ■  s) 

Thermal  conductivity 

J/(m  ■  s  ■  °K) 
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5.10  Input  Variable  Description 

The  variables  that  are  specified  in  the  input  data  file  IDATA  are  defined  below.  These 
variables  are  input  in  several  extended  NAMELIST  blocks.  There  are  defaults  for 
most  input  variables  as  shown  in  the  descriptions  below.  It  is,  therefore,  not  necessary 
to  specify  every  input  variable  in  the  NAMELIST  blocks  in  the  IDATA  file. 


5.10.1  CONTROL  NAMELIST 

TITLE .  A  descriptive  title  consisting  of  up  to  60  characters.  The  character 

string  must  be  delimited  with  single  quotes. 

Internal  Name  =  TITLE 

Default  Value  =  ’DUMPSTER:  AKTEC  DUMP  COMBUSTOR  CODE.’ 

EXECUTION. MODE . . . .  The  mode  of  execution  of  DUMPSTER  may  be  set  to  one  of  the 
following: 

=  ’NEW’  Start  a  new  calculation  from  scratch.  Generate  a 

meoh  and  initialize  the  flow  field. 

=  ’ RESTART’  Start  a  new  calculation  using  the  data  from  a 
restart  file  (IRSTR)  to  initialize  the  flow  field.  The  mesh  in 
the  restart  file  may  aiso  be  used  by  setting  the  input  variable 
MESH. GENERATION. MODE  to  ’RESTART’.  The  flow  field  variables  from 
the  restart  file  will  be  interpolated  on  to  the  mesh  that  is  generated 
for  this  calculation. 

=  ’ CONTINUE’  Continue  a  calculation  saved  on  a  restart  file.  The 
mesh  and  flow  field  variables  from  a  restart  file  (will  be  used  for 
initial  conditions.  The  calculation  will  continue  as  if  the  previous 
calculation  did  not  stop.  The  CFLM,  for  example,  will  begin  at  the 
value  it  had  at  the  end  of  the  previous  calculation  saved  in  the  restart 
file.  The  mesh  cell  indices  (I  and  J)  must  be  the  same  as  those  in 
the  restart  file. 

Internal  Name  =  ICNTRL 
Default  Value  =  ’NEW’ 

NUMBER .  OF .  STEPS .. .  The  number  of  computational  steps  that  the  calculation  will  run. 

If  you  are  restarting  an  old  calculation  with  EXECUTION  .MODE  = 
CONTINUE  then  the  step  counter  (ISTEP)  will  start  at  the  step  in 
the  restart  file,  otherwise  the  step  counter  always  starts  at  0. 

Internal  Name  =  NSTEPS 
Default  Value  =  0 

CONVERGENCE. TOLERANCE  Specify  the  number  of  orders  of  magnitude  that  the  convergence 
level  must  be  reduced  to  achieve  convergence.  The  convergence  level 
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(£„)  of  a  zone  is  defined  as  the  log10  of  the  average  of  the  absolute 
value  of  the  residuals  of  the  mass  density  variable. 

The  average  convergence  level  (CQNVA)  is  simply  the  average  of  the 
zonal  convergence  levels  for  all  of  the  zones. 

c°»v*  =  i  E  e. 

n=l  ,N 

where  N  is  the  total  number  of  zones.  When  CONVA  reduces  by  more 
than  the  CONVERGENCE. TOLERANCE  then  the  calculation  stops. 

Internal  Name  =  CNVTL 
Default  Value  =  -6.0 

CPU .  SECONDS .  MAXIMUM  The  calculation  will  terminate  safely  before  the  number  of  CPU 

seconds  specified  with  this  input  variable  is  reached.  A  restart  file 
will  be  saved.  This  feature  may  not  be  implemented  on  the  version 
of  DUMPSTER  on  your  computer  system. 

Internal  Name  =  SECJQB 
Default  Value  =  999999.0 

NUMBER .  OF .  ZONES . . .  The  number  of  zones  in  the  calculation. 

Internal  Name  =  NZONES 
Default  Value  =  1 

COORDINATE. SYSTEM  Specify  the  type  of  coordinate  system. 

=’ AXISYMMETRIC’  Axisymmetric  (X= axial  and  Y=radius) 

=  ’ PLANAR’  2D  planar 

Internal  Name  =  ICRD 
Default  Value  =  'AXISYMMETRIC' 

DEBUG. FLAGS .  Set  these  inputs  to  turn  on/off  debug  and  standard  out  printout 

statements. 

IDBUG(l)  =  1,  Print  convergence  data  to  standard  out. 

IDBUG(l)  =  2,  Print  convergence  data  to  standard  out  for  each  zone. 
IDBUG(2)  =  1,  Print  pointers  in  subroutine  ARRAYS. 

IDBUG(3)  =  1,  Print  U-array  at  beginning  of  subroutines  FLXRI  and 
FLXRJ. 

Internal  Name  =  (IDBUG(L)  ,L=1 , 10) 

Default  Value  =  0 
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5.10.2  PROPERTIES  NAMELIST 

THERMODYNAMIC. MODEL  Selection  for  the  thermodynamics  model. 

=  ’  PERFECT .  GAS  ’  Thermally  and  calorically  perfect  gas.  You  must 
also  specify  GAMMA  and  GAS. CONSTANT  below. 

=’ THERMALLY . PERFECT . GAS ’  Thermally  perfect  gas  (calorically  im¬ 
perfect)  temperature-dependent  properties.  Also  specify:  GAMMA, 
GAS. CONSTANT,  and  SPECIFIC. HEAT. COEFS 

Internal  Name  =  ITHRMO 
Default  Value  =  ’PERFECT .GAS’ 

GAMMA .  The  ratio  of  specific  heat  capacities  used  in  the  perfect  gas  thermal 

equation  of  state:  P  —  (7  -  1  )pe 

Internal  Name  =  GAMC 
Default  Value  =1.4 

GAS. CONSTANT .  Perfect  gas  constant  as  used  in  the  relation  for  specific  heat  capacity 

Cv  =  R/(i-  1). 

Internal  Name  =  GASC 
Default  Value  =  287.0 

SPECIFIC. HEAT. COEFS  Coefficients  for  temperature-dependent  specific  heat  capacities  at 
constant  volume  for  the  function  shown  below. 

CV(T)  =  C\TC 2  4-  C3T  4-  C4 

the  units  of  the  CVCQN  coefficients  (c,,i  =  1,4)  must  be  choosen  to  be 
consistent  with  the  units  of  temperature  in  degrees- Kelvin  and  the 
units  of  Cv  in  the  units  of  m2/(s 2  •  °/v ). 

Internal  Name  =  (CVCON(L) ,L=1,4) 

Default  Value  =0.0,  1.0,  0.0,  717 . 5  (for  air) 

VISCOSITY. MODEL. . .  Selection  for  the  viscosity  model. 

=  ’INVISCID’  No  viscosity  or  heat  conduction. 

=  ’C0NSTANT’  Viscosity  is  uniform  and  constant.  Also  specify  the 
constant  viscosity  as  the  first  element  in  LAMINAR.  VISCOSITY  .COEFS 

=  ’ LAMINAR’  Temperature-dependent  laminar  viscosity.  Also 
specify  the  viscosity  coefficients  in  LAMINAR. VISCOSITY. COEFS 

= ’TURBULENT’  Use  the  Baldwin/Lomax  algebraic  model  for  eddy 
viscosity  and  temperature-dependent  laminar  viscosity.  Also  specify 
the  laminar  viscosity  coefficients  in  LAMINAR.  VISCOSITY  .COEFS 

=  ’KE. TWO. EQUATION’  k  t  two-equation  model.  Also  specify  the 
k~(  coefficients  in  ’KE. CONSTANTS’  and  the  laminar  coefficients  in 
LAMINAR . VISCOSITY . COEFS 
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Internal  Name  =  IVISC 
Default  Value  =  ’INVISCID’ 

THIN. LAYER. OPTION  The  thin-layer  option  is  an  approximation  to  the  diffusion  terms. 

When  this  option  is  on,  the  cross  derivatives  in  the  diffusion  terms 
are  elimenated. 

=YES  Use  only  the  thin-layer  diffusion  terms. 

=N0  Use  full  diffusion  terms. 

Internal  Name  =  ITHINL 
Default  Value  =  ’NO’ 


LAMINAR. VISCOSITY. COEFS  Coefficients  for  temperature-dependent  laminar  viscosity  (fit) 
are  used  in  the  following  formula. 


W(T)  = 


ClTc*  +  c3T  +  C4 
(c$T  +  c$) 


The  units  of  the  coefficients  (c,,i  =  1,6)  must  be  choosen  to  be  con¬ 
sistent  with  the  units  of  temperature  in  degrees-Kelvin  and  the  units 
of  viscosity  of  I(g/(m  ■  s). 

Internal  Name  =  (VISCON(L) ,L*1 ,6) 

Default  Value  =  1.4519E-6,  1.5,  0.0,  0.0,  1.0,  110.0 

(Sutherland’s  relation  for  air) 

KE. CONSTANTS .  The  constants  used  in  the  k-(  turbulence  model  are  specified  in  this 

array  as  follows. 

TKE(l)  =  C(  1,  Used  in  the  source  term  of  the  c  equation. 

TKE(2)  =  Ct 2,  Used  in  the  source  term  of  the  c  equation. 

TKE(3)  =  Cp,  Used  in  the  definition  of  the  eddy  viscosity. 

TKE(4)  =  Ok,  Used  in  the  diffusion  term  of  the  k  equation. 

TKE(5)  =  (7e,  Used  in  the  diffusion  term  of  the  e  equation. 

Internal  Name  =  (TKECON(L)  ,L*1 ,5) 

Default  Value  =  1.44,  1.92,  0.09,  1.0,  1.3 


CONDUCTIVITY  .MODEL  Selection  for  the  thermal  conductivity  model. 

=  ' CONSTANT’  Enter  constant  as  CONDUCTIVITY. C0EFS(1) 

=  ’PRANDTL. NUMBERS’  Use  laminar  and  turbulent  Prandtl  num¬ 
bers  to  compute  conductivities  from  viscosities.  Also  specify 
LAMINAR. PRANDTL. NUMBER  and  TURBULENT. PRANDTL. NUMBER 
=  ’ TEMPERATURE .  DEPENDENT  ’  Use  turbulent  Prandtl  number  in  cal¬ 
culating  turbulent  conductivity,  and  conductivity  coefficients  to  cal¬ 
culate  the  laminar  conductivity.  Also  specify  CONDUCTIVITY. COEFS 
and  TURBULENT. PRANDTL. NUMBER. 

Internal  Name  =  ICOND 

Default  Value  =  ’PRANDTL. NUMBERS’ 
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LAMINAR. PRANDTL. NUMBER  Laminar  Prandti  number.  May  be  used  to  compute  the  thermal 
conductivity  from  the  viscosity  when  the  thermal  conductivity  model 
(CONDUCTIVITYY. MODEL  is  set  equal  to  ' PRANDTL . NUMBERS *. ) 

Internal  Name  =  PRDL 
Default  Value  =0.71 

TURBULENT . PRANDTL . NUMBER  Turbulent  Prandti  number.  May  be  used  to  compute  the 
turbulent  thermal  conductivity  from  the  viscosity. 

Internal  Name  =  PRDT 
Default  Value  =0.9 

CONDUCTIVITY .  COEFS  The  coefficients  for  temperature-dependent  thermal  conductivity  are 
used  in  the  following  formula. 

.  (T,  _  (ciT*  +  c3T  +  c4) 

(  }  (c5T  +  c6) 

For  Sutherland’s  relation  for  air  as  perfect  gas,  use  the  default  val¬ 
ues  of  the  coefficients  and  also  set  GASC=287.0,  GAMMA=1.4,  and 
PRDL=0.7l.  The  coefficients  must  be  choosen  to  be  consistent  with 
the  units  of  temperature  in  degrees- Kelvin  and  the  units  of  conduc¬ 
tivity  of  Kg  ■  m/°K. 

Internal  Name  =  (CNDCON(L)  ,L=1 ,6) 

Default  Value  =  2.0541E-3,  1.5,  0.0,  0.0,  1.0,  110.0 

NUMBER .  OF . SPECIES  The  number  of  chemical  species  calculated.  This  number  must  not 
be  greater  than  MAXNSPC  (see  Section  5.4).  Trace  chemical  species 
have  no  effect  upon  the  other  species  or  the  fluid  dynamics  of  the 
flow. 

Internal  Name  =  NSPC 
Default  Value  =  0 

SPECIES . COEFS .  Schmidt  numbers  used  to  calculated  the  diffusion  coefficients  for  the 

chemical  species  from  the  viscosity  coefficient  of  the  flow.  The  dif¬ 
fusion  coefficient,  T is  calculated  by  dividing  the  viscosity  Ym  by 
the  Schmidt  number,  a jt,  of  the  -species. 

*7.  = 

Internal  Name  =  (SPCCON(L)  ,L=1  ,NSPC) 

Default  Value  =  0 
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5.10.3  MESH  NAMELIST 


MESH. GENERATION. MODE  The  mesh  generation  mode. 

=  ’ INTERNAL’  Use  the  internal  mesh  gener¬ 

ator.  If  EXECUTION. MODE  =  ’NEW’,  an  internal  mesh  is  generated. 
If  EXECUTION. MODE  =  ’RESTART’  or  ’CONTINUE’,  the  mesh  is  read 
from  the  input  restart  file  IRSTR. 

=  ’ EXTERNAL’  Read  an  external  mesh  from  the  file  IMESH.  DUMP¬ 
STER  checks  the  number  of  cells  in  the  i  and  j  directions  of  the  exter¬ 
nal  mesh  (not  counting  the  boundary  cells).  They  should  be  the  same 
as  the  number  of  cells  read  from  the  IDATA  file  as  NUMBER .  I .  CELLS 
and  NUMBER. J. CELLS. 

Interna]  Name  =  IMSHG 
Default  Value  =  ’INTERNAL’ 
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5.10.4  OUTPUT  NAMELIST 


QUICK  .PRINT. FREQUENCY  The  step  interval  for  printing  the  ‘quick  print  out’.  This  is  a  short 
listing  of  a  the  times  step,  CFLM,  and  a  listing  of  a  few  variables 
along  the  zone  walls.  This  prints  to  the  screen  (Unix  standard  out). 

Internal  Name  =  IFQPX 
Default  Value  =  99999 

NUMBER. OF. PRINTS  .  Number  of  minor  prints  to  make  at  the  computational  steps  defined 
in  STEPS. TO. PRINT.  A  minor  print  lists  various  field  variables  along 
each  zone  boundary  and  some  along  the  middle  row  of  cells  ( J*JD/2). 
Also  the  field  variable  in  zero  or  more  columns  of  cells  (i=constant), 
may  be  printed  out  depending  upon  the  input  values  of  set  in  the 
ZONE. OUTPUT  NAMELIST  for  each  zone. 

Internal  Name  =  NPRT 
Default  Value  =  0 

STEPS.  TO.  PRINT. . . .  This  is  an  array  of  computational  steps  at  which  minor  printouts  are 
done.  There  are  NPRT  values. 

Internal  Name  =  IPRT 
Default  Value  =  0 

NUMBER. OF. PLOTS. . .  Number  of  plots  to  make  at  the  computational  steps  defined  in 
STEPS. TO. PLOT. 

Internal  Name  =  NPLT 
Default  Value  =  0 

STEPS. TO. PLOT .  This  is  an  array  of  computational  steps  at  which  a  plot  is  made. 

These  plots  are  in  TECPLOT  format  and  appended  together  in  the 
output  file  OPLTF.  The  last  plot,  is  done  automatically. 

Internal  Name  =  IPLT 
Default  Value  =  0 

MA  JO!  PRINT .  SKIP .  I  At  the  end  of  a  calculation  all  of  the  field  variables  in  all  of  the  cells  is 
pr.iited  to  OPRNT.  Since  the  amount  of  data  in  a  major  print  out  can 
be  very  large,  used  MAJOR . PRINT. SKIP .  I  and  MA JOR. PRINT. SKIP.  J 
to  skip  points  that  are  printed  out. 

Internal  Name  =  ISKP 
Default  Value  =  999 

MA  JOR.  PRINT.  SKIP.  J  Same  as  MAJOR.  PRINT.  SKIP .  I  except  this  applies  to  the  ;-index  in 
the  major  print  out. 

Internal  Name  =  JSKP 
Default  Value  =  999 


PRINT. PROPERTIES  . 


PRINT. METRICS 


Selection  for  printing  properties  in  the  minor  printouts. 

=  ’N0'  Do  not  print  properties. 

“’YES’  Print  properties. 

Internal  Name  =  IPRTQ 
Default  Value  =  ’NO* 

Selection  for  printing  mesh  coordinates  and  metrics. 

=  ’N0’  Do  not  print  coordinates  and  metrics. 

=  ’XY’  Print  x  and  y  coordinates  only. 

“’VOLUMES'  Print  volume,  and  NAZ  metrics  only. 

=  ’ALL’  Print  i,  y,  volume,  and  all  the  mesh  metics. 

Internal  Name  =  IPRTM 
Default  Value  =  ’NO’ 
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5.10.5  NUMERICS  NAMELIST 


TS.  RELAX  ATI  ON.  FACTOR  Relaxation  coefficient  for  the  overal  time  step.  This  is  a  factor 

between  0  and  1  that  multiplies  DU  before  it  is  added  to  U.  This 
should  not  be  used  for  explicit  calculations.  It  is  useful  for  damping 
implicit  calculations. 

Internal  Name  =  RFTS 
Default  Value  =  1 . 0 

LU.BETA .  This  is  a  coefficient  used  in  the  LU-SGS  relaxation.  It  is  the  ft  factor 

in  Eq.  22.  It  should  normally  be  set  to  1.0.  If  the  LU-SGS  has 
stability  problems  increase  LU.BETA  up  to  a  maximum  of  about  5.0. 

Internal  Name  =  BETA 
Default  Value  =1.0 

TIME .  STEP .  MODE _  This  is  the  method  used  to  calculate  the  base  time  steps.  For  each 

cell  the  base  time  step  is  multiplied  by  the  CFLM  factor  to  get  the 
actual  time  step  for  that  cell.  The  base  time  step  is  calculated  using 
one  of  the  following  options. 

=’ LOCAL’  The  local  explicit  time-step. 

= ’UNIFORM’  A  uniform  minimum  explicit  time  step. 

=’ LOCAL . J . COLUMN ’  A  time  step  is  calculated  for  each  column  of 
cells  (same  value  of  I).  The  time  step  is  explicit  in  the  i -direction.  It 
may  be  larger  than  the  explicit  limit  in  the  j  direction. 

= ’UNIFORM. NEGLECT. J’  This  is  al¬ 

most  the  same  as  ’  LOCAL .  J .  COLUMN  ’  except  that  the  time  step  is 
uniform  throughout  the  zone. 

= 'UNIFORM. NEGLECT.  I’  This  is  the  same  as  ’ UNIFORM. NEGLECT.  J’ 
except  explicit  in  the  j-direction. 

Internal  Name  =  ICDT 
Default  Value  =  ’LOCAL’ 

DU. CHANGE. MAXIMUM  The  relative  changes  in  the  mass  density  (Ap)  and  the  total  internal 
energy  density  (A E)  in  one  time  step  can  be  limited  to  be  no  larger 
than  DU. CHANGE. MAXIMUM.  If  either  exceeds  this  allowable  maximum 
relative  change,  then  the  factor  necessary  to  reduce  them  to  the  limit 
is  also  applied  to  the  other  elements  of  the  U  array.  Normally  this 
option  is  not  used.  Setting  it  to  1.0  turns  it  off.  Use  only  in  cases  of 
severe  numerical  instabilities. 

Internal  Name  =  CHGMAX 
Default  Value  =1.0 

CFLM. BEGIN .  This  is  the  beginning  CFL  multiplier  (CFLM)  for  a  run.  CFLM 

multiplies  the  base  time  step  at  each  cell  (see  the  input  variable 
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TIME .  STEP .  MODE)  to  obtain  the  actual  time  step  to  use  in  the  calcu¬ 
lation.  On  restart  runs  the  CFLM .  BEGIN  is  usually  ignored  and  the  ini¬ 
tial  CFLM  is  set  equal  to  the  value  from  the  restart  file.  If,  however, 
a  minus  sign  is  placed  before  the  value  assigned  to  (CFLM. BEGIN), 
then  the  CFLM  is  set  equal  to  absolute  value  of  CFLM. BEGIN. 

Internal  Name  =  CFLBEG 
Default  Value  =1.0 

CFLM. FACTOR .  This  factor  multiplies  the  CFL  multiplier  (CFLM)  on  each  compu¬ 

tational  step  during  the  calculation  causing  CFLM  to  geometrically 
increase. 

CFLMn  =  CFLM .  FACTOR  *  CFLMn_1 

Internal  Name  =  CFLFAC 
Default  Value  =  1.0 

CFLM. MAXIMUM .  This  is  the  maximum  value  that  the  CFL  multiplier  is  allowed  to  be. 

Internal  Name  =  CFLMAX 
Default  Value  =  1000000.0 

FLUX. FUNCTION. TYPE  DUMPSTER  allows  you  to  select  form  the  following  methods  for 

calculating  the  inviscid  fluxes  at  the  faces  of  computational  cells. 

=  ,R0E’  Roes’s  flux-difference  flux  function. 

=’ STEGER. WARMING’  Steger- Warming  flux  splitting. 

=  ’  HARTEN .  YEE  ’  Harten  and  Yee  flux-difference  flux  function. 

Internal  Name  =  IFLXF 

Default  Value  =  ’STEGER. WARMING’ 

FLUX. FUNCTION. C0EFS  These  coefficients  are  used  in  the  flux  functions  above.  The  defini¬ 
tions  of  the  coefficients  and  their  default  values  depend  upon  which 
flux  function  is  selected.  See  Section  3 


FLUX . FUNCTION . COEFS(L) 

Flux  Function 

1 

2 

3j 

mm 

5 

6 

7 

ROE 

Default  Values 

■ 

■ 

■ 

■ 

Cl 

0.05 

C2 

0.05 

C3 

0.05 

STEGER. WARMING 

Default  Values 

<t>  t 

1.0 

1.0 

C5 

0.05 

C6 

0.05 

C7 

0.05 

HARTEN. YEE 

Default  Values 

■ 

■ 

■ 

■ 

C5 

0.05 

C6 

0.05 

C7 

0.05 

Internal  Name  =  (FLXCFS(L)  ,L*1 ,7) 
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5.10.6  ZONE. INITIAL. CONDITIONS  NAMELIST 


ZONE. NUMBER 


IC. METHOD 


PRESSURE 


TEMPERATURE 


U. VELOCITY 


V. VELOCITY 


W. VELOCITY 


TURBULENT. ENERGY  . 


Number  of  the  zone  for  this  NAMESLIST  block.  The  zones  must  be 
numbered  from  bottom  to  top  starting  with  the  number  1. 

Internal  Name  =  NZ 
Default  Value  =  1 

This  selects  the  method  in  which  the  flow  field  is  initialized  for  this 
zone.  The  flow  field  conditions  are  given  below. 

=  ’UNIFORM . CONDITIONS ’  Uniform  conditions. 

=  ’ ID. NOZZLE*  One-dimensional  nozzle  solution.  The  flow  through 
this  zone  is  calculated  as  if  it  were  isentropic  and  one-dimensional. 
The  PRESSURE  and  TEMPERATURE  are  taken  as  total.  The  Mach 
number  at  the  minimum  flow  area  of  the  channel  is  specified  as 
THRO AT. MACH. NUMBER. 

= ’USER. SUBROUTINE*  User  supplied  subroutine  called  USERIC  is 
supplied  for  the  user  to  add  custom  coding  for  initializing  the  flow 
in  this  zone. 

Internal  Name  =  IIC(NZ) 

Default  Value  =  ’UNIFORM’ 

May  be  static  or  total  pressure  used  for  initializing  the  flow  field  of 
this  zone  depending  upon  the  selection  for  IC. METHOD. 

Internal  Name  =  PRESI(NZ) 

Default  Value  =  101325.0 

May  be  static  or  total  temperature  used  for  initializing  the  flow  field 
of  this  zone  depending  upon  the  selection  for  IC  .METHOD  for  this  zone. 

Internal  Name  =  TEMPI  (NZ) 

Default  Value  =  100.0 

x-component  of  velocity  used  to  initialize  the  (low  field  of  this  zone. 

Internal  Name  =  WELI(l.NZ) 

Default  Value  =  100.0 

r-component  of  velocity  used  to  initialize  the  flow  field  of  this  zone. 

Internal  Name  =  VVELIC2.NZ) 

Default  Value  =0.0 

Swirl-component  of  velocity  used  to  initialize  the  flow  field  of  this 
zone. 

Internal  Name  =  WELI(3,NZ) 

Default  Value  =0.0 

The  turbulent  kinetic  energy  k  used  to  initialize  this  zone  if  the  k-( 
turbulence  model  is  selected. 
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Internal  Nam  3  =  TKI(NZ) 
Default  Value  =0.0 


TURBULENT.  DISSIPATION  The  turbulent  dissipation  rate  e  used  to  initialize  this  zone  if  the 

k-e  turbulence  model  is  selected. 

Internal  Name  =  TEI(NZ) 

Default  Value  =0.0 

MASS. FRACTIONS _  The  mass  fractions  of  the  chemical  species  initialized  in  this  zone. 

Enter  one  value  for  each  species. 

Internal  Name  =  (FMI(N,NZ)  ,N=1,NSPC) 

Default  Value  =0.0 

THROAT .  MACH .  NUMBER  When  the  ONE .  DIMENSIONAL .  NOZZLE  is  selected  for  the  initialization 

method,  this  variable  specifies  the  Mach  number  at  the  minimum 
area  of  the  zone  assuming  flow  in  the  i-direction. 

Internal  Name  =  THMACH(NZ) 

Default  Value  =  1.0 
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5.10.7  ZONE. BOUNDARY. CONDITIONS  NAMELIST 


In  this  NAMELIST  block  the  boundary  conditions  are  selected  for  the  six  segments 
of  a  zone.  Some  boundary  condition  types  may  only  be  selected  for  certain  boundary 
segments.  The  allowable  boundary  condition  types  for  each  boundary  segment  are 
tabulated  below. 


Zone  Boundaries 

Also  Specify: 

B.C.  Type 

Top 

Bottom 

Left 

Right 

'FREE. SLIP. WALL’ 

mm 

mm 

Y 

mm 

’NO. SLIP. WALL’ 

1 

Y 

I 

’WALL. FUNCTION’ 

I 

1 

Y 

O 

’INTERZONE’ 

H 

■1: 

N 

N 

’SUBSONIC. INFLOW’ 

N 

N 

Y 

N 

’SUPERSONIC. INFLOW’ 

N 

N 

Y 

N 

Static  pres¬ 

sure,  static  temperature,  and 
velocity  components 

’SUBSONIC. OUTFLOW’ 

N 

N 

N 

Ambient  static  pressure  (PEX) 

’ SUPERSONIC . OUTFLOW ’ 

N 

N 

N 

’MIXED. OUTFLOW’ 

N 

N 

N 

Ambient  static  pressure 

’PERIODIC’ 

N 

N 

Y 

’VORTEX. OUTFLOW’ 

N 

N 

N 

Ambient  static  pressure 

For  supersonic  inflows  specify  static  conditions  for  all  variables.  For  subsonic  inflows 
specify  total  conditions  for  pressure  and  temperature  and  specify  the  respective  di¬ 
rection  cosines  for  the  velocities  such  that  the  square  root  of  the  sum  of  the  squares 
of  the  direction  cosines  is  equal  to  1.0. 


ZONE. NUMBER .  Number  of  the  zone  for  this  NAMESLIST  block.  The  zones  must  be 

numbered  from  bottom  to  top  starting  with  the  number  1. 

Internal  Name  =  NZ 
Default  Value  =  1 


BC. TYPE. BOTTOM. SI  Boundary  condition  type  on  the  first  segment  on  the  bottom  bound¬ 
ary  of  the  zone. 

Internal  Name  =  KBCB(l.NZ) 

Default  Value  =  ’FREE. SLIP. WALL’ 

BC. TYPE. BOTTOM. S2  Boundary  condition  type  on  the  second  segment  on  the  bottom 
boundary  of  the  zone. 

Internal  Name  =  KBCB(2,NZ) 

Default  Value  =  'FREE. SLIP. WALL’ 


BC  .TYPE  .BOTTOM  .S3  Boundary  condition  type  on  the  third  segment  on  the  bottom  bound¬ 
ary  of  the  zone. 
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Usually 

Inflow 


KBCT(1) 


T" 

IBCT(2) 


KBCT(2) 


T" 

IBCT  (3 ) 


KBCT(3) 


KBCL 


Lr- Zone  Boundary 


Zone  NZ 


Usually 

Oulllow 

KBCR 


IBCB<2)  IBCB(3) 

l _ ,  1, 

—  l  III  . I  ill 

KBCB(1)  - — KE>OB(2>  - -  KBCB(3) 


Figure  8:  Zone  Boundary  Condition  Segments  Shown  in  Computational  Space. 


Internal  Name  =  KBCB(3,NZ) 

Default  Value  =  FREE .  SLIP .  WALL  * 

BC. TYPE. TOP. SI ....  Boundary  condition  type  on  the  first  segment  on  the  top  boudaxy  of 

the  zone. 

Internal  Name  =  KBCT(1,NZ) 

Default  Value  =  ’FREE. SLIP. HALL’ 


BC  .TYPE .  TOP .  S2 _  Boundary  condition  type  on  the  second  segment  on  the  top  boundary 

of  the  zone. 

Internal  Name  =  KBCT(2,NZ) 

Default  Value  =  'FREE. SLIP. WALL’ 

BC. TYPE. TOP. S3 Boundary  condition  type  on  the  third  segment  on  the  top  boundary 

of  the  zone. 

Internal  Name  =  KBCT(3,NZ) 

Default  Value  =  ’FREE. SLIP  .WALL’ 

BC. TYPE. LEFT .  Boundary  condition  type  on  the  left  boundary  of  the  zone.  This  is 

usually  the  inflow  boundary. 

Interna]  Name  =  KBCL(NZ) 

Default  Value  =  ’SUBSONIC. INFLOW’ 
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BC . TYPE .  RIGHT .  Boundary  condition  type  on  the  right  boundary  of  the  zone.  This  is 

usually  the  outflow  boundary. 

Internal  Name  =  KBCR(NZ) 

Default  Value  =  ’SUPERSONIC . OUTFLOW* 

AMBIENT. PRESSURE  .  This  is  the  static  pressure  used  for  setting  the  outflow  boundary 
conditions. 

Internal  Name  =  PEX(NZ) 

Default  Value  =  101325.0 

AMBIENT . TEMPERATURE  This  is  the  static  temperature  that  is  used  for  setting  outflow 

boundary  conditions.  This  is  only  used  for  the  case  when  the  flow 
reverses  and  flows  into  the  boundary. 

Internal  Na^e  =  TEX(NZ) 

Default  Value  =  273.15 

WALL . TEMPERATURE . LEFT  This  is  the  wall  temperature  on  the  left  boundary  used  for  com¬ 
puting  the  surface  heat  transfer  rate. 

=0.0  If  set  to  zero,  the  boundary  is  adiabatic. 

=greater  thanO .  0  If  nonadiabatic,  set  equal  to  the  wall  temperature. 

Internal  Name  =  TWL(NZ) 

Default  Value  =0.0 

WALL. TEMPERATURE. RIGHT  Same  as  above  for  the  rig!:'  boundary. 

Internal  Name  =  TWR(NZ) 

Default  Value  =0.0 

WALL . TEMPERATURE .  BOTTOM  Same  as  above  for  the  bottom  boundary. 

Internal  Name  =  TWB(NZ) 

Default  Value  =0.0 

WALL  .TEMPERATURE. TOP  Same  as  above  for  the  top  boundary. 

Internal  Name  =  TWT(NZ) 

Default  Value  =0.0 

BC.I .  INDEX  .BOTTOM  The  bottom  boundary  of  the  zone  is  divided  into  three  segments  as 
shown  in  Figure  8.  A  different  boundary  condition  type  may  be  spec¬ 
ified  for  each  boundary  segme  it.  The  boundary  segments  are  defined 
by  the  ending  i-index  for  each  segment.  BC.I  .INDEX. BOTTOM (1)  is 
the  ending  i-index  for  segment  one  and  BC .  I .  INDEX  .BOTTOM (2)  is 
the  ending  i  index  for  segment  two. 

Internal  Name  =  (IBCB(L ,NZ)  ,L=2 ,3) 

Default  Value  =2,2 


BC .  I .  INDEX .  TOP _  Same  as  above  for  the  top  boundary. 

Internal  Name  =  (IBCT(L.NZ)  ,L*2,3) 

Default  Value  =  2,2 

INTERPOLATION  .METHOD  The  inflow  boundary  conditions  on  the  left  boundary  may  be  spec¬ 
ified  by  a  profile  for  each  variable.  This  profile  is  defined  in  a  table  of 
data  points.  The  boundary  conditions  at  a  cell  are  calculated  from 
these  data  points  by  interpolation.  The  method  used  to  interpolate 
from  the  input  data  table  to  the  boundary  cells  may  be  linear  or 
quadratic  as  specified  here. 

=  ’ LINEAR’  Linear  interpolation  between  data  points. 

=  ’ QUADRATIC’  Quadratic  interpolation  between  data  points. 
Quadratic  interpolation  can  result  in  interpolated  points  being  much 
different  from  the  values  of  the  source  data  points  from  which  they 
are  interpolated.  This  is  especially  when  the  source  data  points  are 
not  smooth.  When  transitioning  from  one  smooth  region  to  another 
smooth  region  use  several  points  very  close  together. 

Internal  Name  =  NTRPBC(NZ) 

Default  Value  =  ’LINEAR* 

NPTS.UVWPT .  ARRAY  .  Number  of  data  points  defining  the  inflow  boundary  condition  profile 
on  left  boundary.  If  only  one  data  point  is  defined,  then  a  uniform 
profile  is  assumed. 

Internal  Name  =  NBCVPT(NZ) 

Default  Value  =  1 

UVWPT .  ARRAY .  This  is  an  array  that  contains  data  points  defining  an  inflow  boundary 

condition  profile  at  the  left  boundary  of  the  zone.  The  internal  name 
for  this  array  is  BCVPT(L,N,NZ)  It  has  three  subscripts;  the  first  is 
for  the  data  point  number  (L),  the  second  indicates  the  variable  (N) 
which  varies  from  1  to  NBCVPT(NZ),  and  the  third  indicates  the  zone 
number  (NZ). 

BCVPT(L ,  1  ,NZ)  =  y-  or  r-coordinate 
BCVPT(L,2,NZ)  =  x- velocity 
BCVPT(L  ,3  ,NZ)  =  y- velocity 
BCVPT(L,4,NZ)  =  z-velocity 
BCVPTCL ,5 ,NZ)  =  pressure 
BCVPT(L,6,NZ)  =  temperature 

The  default  values  are  zero  for  all  elements. 

Internal  Name  =  (BCVPT(L,N,NZ) ,N=1 ,6) ,L*1 .NBCVPT(NZ) 
Default  Value  =0.0 


NPTS.KE. ARRAY 


Number  of  data  points  defining  the  k-e  turbulence  variables  for  the 
inflow  boundary  condition  profile  on  the  left  boundary. 


Internal  Name  =  NBCTKE(NZ) 

Default  Value  =  1 

KE. ARRAY .  This  array  contains  the  data  points  for  defining  the  k-(  turbulence 

variables  for  the  inflow  boundary  condition  on  the  left  boundary  of 
zone  NZ.  The  internal  name  for  this  input  variable  is  BCTKE(L,N,NZ) . 
The  first  subscript  denotes  the  data  point  number  (ranging  from  1 
to  NPTS.KE. ARRAY),  the  second  indicates  the  variable,  and  the  third 
indicates  the  zone  number. 

BCTKE(L,1,NZ)  =  y-coordinate  ( y ) 

BCTKE(L,2,NZ)  =  turbulence  energy  ( k ) 

BCTKE(L,3,NZ)  =  turbulence  dissipation  (c) 

The  default  values  are  zero  for  all  elements. 

Internal  Name  =  (BCTKE(L,N,NZ)  ,N=1,2)  ,L=1 .NBCTKE(NZ) )  ‘ 
Default  Value  =  0.0 

NPTS.S. ARRAY .  Number  of  data  points  defining  the  chemical  species  profile  for  the 

inflow  boundary  condition  on  the  left  boundary. 

Internal  Name  =  NBCSPC(NZ) 

Default  Value  =  1 

S. ARRAY .  This  array  contains  the  data  points  for  defining  the  chemical  species 

profile  for  the  inflow  boundary  condition  on  the  boundary  of  zone 
NZ.  The  internal  name  for  this  input  variable  is  BCSPC(L.N.NZ).  The 
first  subscript  denotes  the  data  point  number  (ranging  from  1  to 
NPTS.S. ARRAY),  the  second  indicates  the  variable,  and  the  third  in¬ 
dicates  the  zone  number. 

BCSPCfL,  1  ,NZ)  =  y-coordinate 
BCSPC(L,2,NZ)  =  Mass  fraction  for  specie  1  (/]) 
BCSPC(L,3,NZ)  =  Mass  fraction  for  specie  2  (/2) 

BCSPC(L,NSPC-*T  ,NZ)  =  Mass  fraction  for  species  NSPC 

The  default  values  are  zero. 

Internal  Name  = 

((BCSPC(L.N.NZ) ,N=1 .NSPC+l) ,L=1 ,NBCSPC(NZ) ) 

Default  Value  =0.0 
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5.10.8  ZONE. GEOMETRY  NAMELIST 


The  geometry  definition  arrays  may  extend  beyond  the  actual  boundaries  of  the  mesh 
that  is  to  be  generated.  The  geometry  arrays  must  be  single- valued  functions  of  x. 
For  nozzle  and  dump  combustor  calculations,  the  mesh  must  be  oriented  so  that  the 
flow  is  directed  in  the  positive  x-direction. 

If  NUMBER. OF  .POINTS  .BOTTOM  =  0  then  the  geometry  of  the  bottom  surface  is  ob¬ 
tained  from  the  user  supplied  subprogram  function  YGEOM(NZ,IBT,XN),  where:  NZ  is 
the  zone  number,  IBT  is  the  integer  control  specifying  bottom  (IBT=l)  or  top  (IBT=2) 
boundary,  and  XN  is  the  x-coordinate  that  is  input  to  the  function.  The  r-coordinate 
(or  y-coordinate)  value  is  returned  as  the  value  of  YGEOM. 


ZONE.  NUMBER .  Number  of  the  zone  for  this  NAMESLIST  block.  The  zones  must  be 

numbered  from  bottom  to  top  starting  with  the  number  1. 

Internal  Name  =  NZ 
Default  Value  =  1 

CONVERSION. FACTOR  The  input  coordinates  for  wall  geometry  are  multiplied  by 

CONVERSION. FACTOR  before  being  used.  This  input  variable  can  be 
used  for  converting  the  units  of  the  input  geometry  data  to  units  of 
meters. 

Internal  Name  =  XYCONV(NZ) 

Default  Value  =1.0 

NUMBER. OF. POINTS. BOTTOM  The  number  of  data  points  defining  the  geometry  of  the  bot¬ 
tom  of  the  zone. 

Internal  Name  =  NSB(NZ) 

Default  Value  =  2 

INTERPOLATION. METHOD. BOTTOM  Selection  of  the  method  used  for  interpolating  the  actual 
geometry  from  the  geometry  data  points  input  with  the  X. BOTTOM 
and  Y. BOTTOM  arrays.  See  notes  on  linear  and  quadratic  in¬ 
terpolation  for  the  INTERPOLATION  .METHOD  input  variable  in  the 
ZONE. BOUNDARY. CONDITIONS  NAMELIST  block. 

=  ’ LINEAR’ 

=’ QUADRATIC’ 

Internal  Name  =  NTRPGB(NZ) 

Default  Value  =  ’LINEAR’ 

X. ORIGIN. BOTTOM. . .  This  is  the  origin  of  the  coordinate  system  used  for  inputting  the 
geometry  for  the  bottom  of  the  zone.  The  final  x  and  y  coordinates 
used  in  the  code  are  equal  to  the  sum  of  the  input  values  x  and  y 
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plus  the  x  and  y  origins  ( xa  and  y0)  as  follows. 


x  -  x0  +  x 

y  =  Vo  +  y 

Internal  Name  =  XSTO(NZ) 

Default  Value  =0.0 

Y. ORIGIN. BOTTOM. . .  The  y  origin  for  the  geometry  of  the  bottom  of  the  zone. 

Internal  Name  =  YSTO(NZ) 

Default  Value  =0.0 

X.  BOTTOM .  This  is  a  vector  of  values  that  define  the  x-coordinates  (i)  of  the 

bottom  geometry  of  the  2one.  The  actual  geometry  that  is  used  is 
interpolated  form  these  values  and  may  be  shifted  and  scaled  using 
CONVERSION. FACTOR  and  X. ORIGIN. 

Internal  Name  =  (XSB(L ,NZ) ,L=1 ,NSB(NZ) ) 

Default  Value  =  0.0, 2.0 

Y.  BOTTOM .  This  is  a  vector  of  values  that  define  the  y  (or  r)  coordinates  ( y )  of 

the  bottom  geometry  of  the  zone.  The  actual  geometry  that  is  used 
is  interpolated  form  these  values  and  may  be  shifted  and  scaled  using 
CONVERSION. FACTOR  and  Y. ORIGIN. 

Internal  Name  =  (YSB(L.NZ)  ,L=1  ,NSB(NZ) ) 

Default  Value  =  0.0, 2.0 

NUMBER. OF. POINTS. TOP  Same  as  NUMBER. OF. POINTS. BOTTOM  except  for  the  geometry  at 
the  top  of  the  zone. 

Internal  Name  =  NST(NZ) 

Default  Value  =  2 

INTERPOLATION .  METHOD .  TOP  Same  as  NUMBER .  OF .  POINTS .  BOTTOM  except  for  the  geometry 
at  the  top  of  the  zone. 

=’ LINEAR’ 

=’ QUADRATIC’ 

Internal  Name  =  NTRPGT(NZ) 

Default  Value  =  'LINEAR' 

X.  ORIGIN. TOP .  Same  as  NUMBER .  OF  .POINTS . BOTTOM  except  for  the  geometry  at  the 

top  of  the  zone. 

Internal  Name  =  XSTD(NZ) 

Default  Value  =0.0 

Y.  ORIGIN  .TOP .  Same  as  NUMBER .  OF  .POINTS. BOTTOM  except  for  the  geometry  at  the 

top  of  the  zone. 
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Internal  Name  =  YSTO(NZ) 

Default  Value  =0.0 

X .  TOP .  Same  as  NUMBER .  OF . POINTS. BOTTOM  except  for  the  geometry  at  the 

top  of  the  zone. 

Internal  Name  =  (XST(L,NZ) ,L*1 ,NSB(NZ) ) 

Default  Value  =  0.0 ,2.0 

Y. TOP .  Same  as  NUMBER. OF  .POINTS. BOTTOM  except  for  the  geometry  at  the 

top  of  the  zone. 

Internal  Name  =  (YST(L.NZ) ,L=1,NSB(NZ)) 

Default  Value  =  0.0, 2.0 
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5.10.9  ZONE. MESH  NAMELIST 

ZONE. NUMBER .  Number  of  the  zone  for  this  NAMESLIST  block.  The  zones  must  be 

numbered  from  bottom  to  top  starting  with  the  number  1. 

Internal  Name  =  NZ 
Default  Value  =  1 

NUMBER .  OF .  CELLS .  I  The  number  of  internal  cells  (not  counting  boundary  cells)  for  this 
zone  in  the  i-direction.  On  the  Cray  computer  and  other  computers 
with  banked  memory  you  should  use  odd  number  of  cells  to  prevent 
bank  conflicts  when  accessing  memory.  This  variable  must  be  consis¬ 
tent  with  the  mesh  generation  input  data  defined  below,  or  consistent 
with  an  externally  generated  mesh. 

Internal  Name  =  IDIM(NZ) 

Default  Value  =  21 

NUMBER . OF . CELLS .  J  The  number  of  internal  cells  (not  counting  boundary  cells)  in  the 
j-direction.  See  discussion  above. 

Internal  Name  =  JDIM(NZ) 

Default  Value  =21 


If  an  external  computational  mesh  is  read  from  the  IMESH  file  or  from  the  restart  file 
IRSTR  (i.e.  MESH. GENERATION. MODE  =  ' EXTERNAL '  or  'RESTART'),  then  the  mesh 
generation  data  below  will  be  ignored. 

The  internal  mesh  generation  routine  creates  meshes  that  have  straight  vertical  i- 
lines  (constant  i  index).  The  mesh  is  generated  in  four  segments  in  the  i-direction 
(which  is  often  the  same  as  the  r-direction)  and  one  segment  in  the  j-direction  a 
schematic  diagram  is  shown  in  Figure  9. 

Do  not  include  boundary  cells  when  calculating  the  following  input  variables. 

The  i-index  at  the  first  interior  cell  on  the  left  boundary  is  three,  since  their  are  two 
boundary  cells.  The  i-index  of  the  interior  cell  that  is  farthest  to  the  right  on  the 
right  boundary  has  an  i-index  of  ID-2.  ID  is  the  total  number  of  interior  and  exterior 
cells  in  the  i-directions.  JD  is  the  total  number  of  interior  and  exterior  cells  in  the 
^-directions.  Subtract  4  from  ID  or  JD  to  calculated  the  number  of  internal  cells  in 
the  i  or  j  direction. 


X. CENTER .  z-coordinate  of  the  line  between  left  and  right  center  regions  as 

shown  in  Figure  9. 
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Internal  Name  =  XCTR(NZ) 

Default  Value  =  0.0 

DELTA. X .  i -direction  mesh  increment  at  X. CENTER. 

Internal  Name  =  DELX(NZ) 

Default  Value  =1.0 

X .  STRETCH .  MODE _  Selection  of  mesh  stretching  mode  in  the  i-direction. 

■  ’  STRETCH .  FACTORS '  Specify  stretch  factors 

for  the  STRETCH . LENGTH . x  (where  x  is  LEFT,  RIGHT,  LEFTCENTER, 
or  RIGHTCENTER)  input  variables  below.  The  mesh  increments  will 
be  stretched  away  from  the  X .  CENTER. 

= ’ LENGTHS ’  Specify  lengths  of  each  region  for  the  STRETCH .  LENGTH 
input  variables  below.  The  mesh  spacing  is  calculated  such  that  the 
number  of  cells  specified  in  the  NUMBER .  OF .  CELLS .  x  (where  x  is  LEFT, 
RIGHT,  LEFTCENTER,  or  RIGHTCENTER)  input  variables  will  extend  the 
specified  length. 

Internal  Name  =  ISW(NZ) 

Default  Value  =  ’LENGTHS' 

NUMBER .  OF .  CELLS .  LEFT  The  number  of  cells  in  the  segment  that  is  on  the  left  boundary. 

Internal  Name  =  NCX(1,NZ) 

Default  Value  =  0 

NUMBER. OF. CELLS. LEFTCENTER  Z  number  of  cells  in  the  segment  that  is  just  left  of 

X.  CENTER. 


Usually 

Outflow 


Figure  9:  Mesh  Generation  Schematic  for  the  x-Direction. 
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Internal  Name  =  NCX(2,NZ) 

Default  Value  =  DEF=0 

NUMBER .  OF .  CELLS .  RIGHTCENTER  The  number  of  cells  in  the  segment  that  is  just  right  of 

X. CENTER. 

Internal  Name  =  NCX(3,NZ) 

Default  Value  =  DEF=0 

NUMBER .  OF .  CELLS .  RIGHT  The  number  of  cells  in  the  segment  that  is  on  the  right  boundary. 

Internal  Name  =  NCX(4,NZ) 

Default  Value  =  0 

STRETCH .  LENGTH .  LEFT  This  is  either  the  stretch  factor  (E)  or  the  length  of  the  LEFT  seg¬ 
ment  depending  upon  X .  STRETCH .  MODE. 

Internal  Name  =  SLX(l,NZ) 

Default  Value  =1.0 

STRETCH . LENGTH . LEFTCENTER  This  is  either  the  stretch  factor  (E)  or  the  length  of  the 
LEFTCENTER  segment  depending  upon  X . STRETCH . MODE. 

Internal  Name  =  SLX(2,NZ) 

Default  Value  =1,0 

STRETCH. LENGTH. RIGHTCENTER  This  is  either  the  stretch  factor  (E)  or  the  length  of  the 
RIGHTCENTER  segment  depending  upon  X. STRETCH. MODE. 

Internal  Name  =  SLX(3,NZ 
Default  Value  =1.0 

STRETCH. LENGTH  .RIGHT  This  is  either  the  stretch  factor  ( E )  or  the  length  of  the  RIGHT 
segment  depending  upon  X. STRETCH. MODE. 

Internal  Name  =  SLX(4,NZ) 

Default  Value  =1.0 

J .  STRETCH .  MODE . . . .  This  is  the  selection  for  the  mode  in  which  the  mesh  cell  spacing  in 
the  vertical  direction  is  calculated. 

*’ STRETCH. FACTORS’  Specify  stretch  factors  and  the  number  of 
blocks  of  j-cells. 

»’ DELY. BOTTOM’  Specify  DELY. BOTTOM  at  the  bottom  of  the  zone 
( j  =  3)  and  the  number  of  j-cells 

in  J. BLOCK. NUMBER. OF. CELLS (1).  The  code  will  then  calculate  a 
stretch  factor  to  fill  the  distance  from  the  bottom  to  the  top  of  the 
zone. 

■’ DELY. TOP’  Specify  DELY. TOP  at  top  of  the  zone  (j  =  JDM)  and 
the  number  of  j-cells  in  J. BLOCK. NUMBER. OF. CELLS (1).  The  code 
will  then  calculate  a  stretch  factor  to  fill  the  distance  from  the  top 
to  the  bottom  of  the  zone. 


=  *  DELY .  TOP .  AND .  BOTTOM  *  Specify  DELY .  TOP  at  top  of  the  zone 
(j  =  JDM),  DELY. BOTTOM  at  the  bottom  of  the  zone,  and  the  num¬ 
ber  of .? -cells  in  J. BLOCK. NUMBER. OF. CELLS (1).  The  code  will  then 
calculate  a  stretch  factor  to  fill  the  distance  between  the  top  and  the 
bottom  of  the  zone. 

Internal  Name  =  JSW(NZ) 

Default  Value  =  *  STRETCH .  FACTORS » 

NUMBER .  OF .  J .  BLOCKS  The  number  of  j-blocks  of  cells  in  zone.  The  maximum  number 

is  5.  The  stretch  factor  is  constant  within  each  ji-block.  The  y- 
increment  in  the  lowest  cell  of  blocks  is  equal  to  the  y-increment  of 
the  cell  below  it  times  the  stretch  factor  in  the  block.  This  variable 
is  ignored  if  J . STRETCH . MODE  is  not  set  to  STRETCH. FACTORS. 

Internal  Name  =  NJBK(NZ) 

Default  Value  =  1 

J .  BLOCK .  NUMBER .  OF .  CELLS  The  number  of  j-cells  in  j-block  L,  where  L  goes  from  1  to 

NUMBER. OF.  J. BLOCKS.  The  first  block  is  at  the  bottom  of  the  zone, 
j-blocks  2  through  NUMBER.  OF .  J  .BLOCKS  are  stacked  above  it. 

Internal  Name  =  (JBK(L.NZ)  ,L=1  ,NJBK) 

Default  Value  =  21 

J. BLOCK. NUMBER. OF. I. COLUMNS  The  number  of  i-ceUs  at  which  the  following  y 

mesh  parameters  are  defined.  The  number  of  jf-blocks 
(NUMBER. OF.  J. BLOCKS)  and  the  number  of  y'-cells  in  each  y-block 
are  the  same  at  each  i-cell.  The  maximum  is  five. 

Internal  Name  =  NXJBK(NZ) 

Default  Value  =  1 

J  .BLOCK .  I  .COLUMNS  The  i-cells  at  which  the  following  j-mesh  parameters  are  defined. 

The  first  z-cell  must  equal  1  and  the  last  must  ID. 

Internal  Name  =  (IJBK(L.NZ)  ,L=1  .NIJBK) 

Default  Value  =  1 

J .  BLOCK .  DELY .  BOTTOM  The  y— cell  increment  at  the  bottom  of  the  zone  that  is  used  to 

generate  the  y-mesh.  This  is  used  if  J . STRETCH . MODE  is  set  to 
DELY. BOTTOM  or  DELY. TOP. AND. BOTTOM. 

Internal  Name  =  (DYBJBK(L.NZ)  ,L=1 , NIJBK) 

Default  Value  =0.0 

J. BLOCK. DELY. TOP  .  The  y  cell  increment  at  the  top  of  zone  that  is  used  to  generate 

the  y-mesh.  This  is  used  if  J . STRETCH . MODE  is  set  to  DELY. TOP  or 
DELY. TOP. AND. BOTTOM. 

Internal  Name  =  (DYTJBK(L.NZ)  ,L=1  .NIJBK) 

Default  Value  =0.0 
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J .  BLOCK .  STRETCH .  FACTORS  Stretch  factors  defined  for  j-blocks  at  the  t-cells  defined  by 

J. BLOCK. I. COLUMNS.  This  variable  is  a  two-dimensional  array.  The 
first  subscript  varies  from  1  to  the  number  of  NI JBK  and  the  second 
subscript  varies  from  1  to  NJBK.  The  stretch  factor  in  each  block  of 
j-cells  is  constant  throughout  the  block.  This  variable  is  ignored  if 
J.  STRETCH.  MODE  is  not  equal  to  STRETCH  .FACTORS. 

Internal  Name  =  ((EJBK(L,  J  ,NZ)  ,L=1  ,NI JBK)  ,  J*1  ,NJBK) 
Default  Value  =  1 
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5.10.10  ZONE. OUTPUT  NAMELIST 


ZONE. NUMBER .  Number  of  the  zone  for  this  NAMESLIST  block.  The  zones  must  be 

numbered  from  bottom  to  top  starting  with  the  number  1. 

Internal  Name  =  NZ 
Default  Value  =  1 

NUMBER. OF.  I . COLUMNS  The  minor  printout  occurs  during  the  calculation  for  a  select  num¬ 
ber  of  columns  of  i-cells.  This  variable  specifies  the  number  of  i-cell 
columns  to  print. 

Internal  Name  =  NCOLPRT(NZ) 

Default  Value  =  0 

I .  COLUMNS .  This  is  an  array  of  i-cell  columns  that  will  be  printed  in  the  minor 

printout. 

Internal  Name  =  (ICOLPRT(L.NZ)  ,L=1  .NCOLPRT(NZ)) 

Default  Value  =  2 

NUMBER. OF. TAPS. . . .  You  may  select  computational  cells  in  the  zone  where  the  flow  field 

values  will  be  monitored  during  the  calculation.  At  the  end  of  the 
calculation  the  tap  values  are  formated  for  input  into  TECPLOT. 
Tap  data  allows  you  to  plot  time  histories  of  various  flow  field  vari¬ 
ables  at  specific  places  in  the  flow  field. 

Internal  Name  =  NTAPS(NZ) 

Default  Value  =  0 

TAP.  I .  INDEX .  This  is  an  array  of  i-indexes  of  the  taps  for  this  zone. 

Internal  Name  =  (ITAP(L.NZ)  ,L=1  ,NTAPS) 

Default  Value  =  1 

TAP .  J  .  INDEX .  This  is  an  array  of  j-indexes  of  the  taps  for  this  zone. 

Internal  Name  =  (JTAP(L,NZ)  ,L=1  ,NTAPS) 

Default  Value  =  1 
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5.'  1  Summary  of  Input  Variables 

The  following  is  a  listing  of  the  input  data  showing  the  format  of  the  input  data.  The 
values  shown  are  the  default  values. 


SCONTROL 

EXECUTION. MODE 

NUMBER. OF. STEPS 

CONVERGENCE . TOLERANCE 

CPU .  SECONDS .  MAXIMUM 

NUMBER. OF. ZONES 

COORDINATE. SYSTEM 

DEBUG. FLAGS 

SEND 

=  "NEW", 

■  o, 

*  -6.0, 

*  99999.9, 

*  1. 

*  "AXISYMMETRIC" , 

=  10*0, 

SPROPERTIES 

THERMODYNAMIC .MODEL 

=  "PERFECT. GAS", 

GAMMA 

=  1.4, 

PERFECT . GAS . CONSTANT 

=  287.0, 

SPECIFIC. HEAT. COEFS 

«  0.0,  1.0,  0.0,  717.5, 

VISCOSITY. MODEL 

*  "INVISCID". 

THIN. LAYER. OPTION 

*  -OFF", 

LAMINAR .VISCOSITY . COEFS 

*  1.4519E-06,  1.5,  0.0, 

0.0,  1.0,  110.0, 

KE. CONSTANTS 

=  1.44,  1.92,  0.09,  1.0 

,  1.3, 

CONDUCTIVITY. MODEL 

*  "NONCONDUCTING", 

LAMINAR . PRANDTL . NUMBER 

=  0.7, 

TURBULENT . PRANDTL . NUMBER 

=  0.9, 

CONDUCTIVITY. COEFS 

=  0.0020541,  1.5,  0.0, 

0.0,  1.0,  110.0, 

NUMBER. OF. SPECIES 

»  0 

SPECIES. COEFS 

=  5*0.0, 

SEND 

SMESH 

MESH . GENERATION . MODE 

=  "INTERNAL", 

ADAPT. STEP. BEGIN 

=  99999, 

ADAPT. STEP. END 

*  99999, 

ADAPT. STEP. INTERVAL 

■  99999, 

ADAPT . RELAXATION . FACTOR 

=  0.2, 

SEND 

$ OUTPUT 

QUICK. PRINT. FREQUENCY  =  99999 

NUMBER. OF. PRINTS  *  0 

STEPS. TO. PRINT  ■  100*0, 
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NUMBER. OF. PLOTS 

STEPS. TO. PLOT 

MAJOR. PRINT. SKIP. I 

MAJOR. PRINT. SKIP. J 

PRINT. PROPERTIES 

PRINT. METRICS 
$END 

*  0 

a  100*0, 

=  999 

=  999 

=  "HO" 

.  "NO" 

$ NUMERICS 

TS . RELAXATION . FACTOR 

=  1.0, 

LU.BETA 

=  1.0, 

TIME. STEP. MODE 

=  "LOCAL", 

DU. CHANGE. MAXIMUM 

=  1.0, 

CFLM. BEGIN 

=  1.0, 

CFLM. FACTOR 

=  1.0, 

CFLM. MAXIMUM 

=  1.0E+09, 

FLUX. FUNCTION. TYPE 

=  "ROE", 

FLUX . FUNCTION .COEFS 

=  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 

$END 

$ZONE . INITIAL . CONDITIONS 

ZONE. NUMBER 

=  1, 

IC. METHOD 

=  "UNIFORM. CONDITIONS", 

RRESSURE 

*  1.01325E+05, 

TEMPERATURE 

*  100.0, 

U. VELOCITY 

=  100.0, 

V. VELOCITY 

=  0.0, 

H. VELOCITY 

=  0.0, 

THROAT. MACH. NUMBER 

=  1.0, 

TURBULENT. ENERGY 

=  0.0, 

TURBULENT . DISSIPATION 

=  0.0, 

MASS. FRACTIONS 

=  0.0, 

$END 

IZONE . BOUNDARY . CONDITIONS 

ZONE. NUMBER 

=  1, 

BC. TYPE. BOTTOM. SI 

=  "FREE. SLIP. HALL", 

BC. TYPE. BOTTOM. S2 

=  "FREE. SLIP. HALL", 

BC. TYPE. BOTTOM. S3 

=  "FREE. SLIP. HALL", 

BC. TYPE. TOP. SI 

=  "FREE. SLIP. HALL", 

BC. TYPE. TOP. S2 

=  "FREE. SLIP. HALL", 

BC. TYPE. TOP. S3 

=  "FREE. SLIP. HALL", 

BC. TYPE. LEFT 

=  "SUBSONIC. INFLOH", 

BC. TYPE. RIGHT 

=  "SUPERSONIC. OUTFLOH", 

AMBIENT. PRESSURE 

=  1.01325E+05, 

AMBIENT . TEMPERATURE 

=  273.15, 

87 


BOUNDARY . CONDITION . COEFS 
WALL . TEMPERATURE . LEFT 
WALL . TEMPERATURE . RIGHT 
WALL . TEMPERATURE . BOTTOM 
WALL . TEMPERATURE . TOP 

BC. I. INDEX. BOTTOM 

BC. I. INDEX. TOP 
INTERPOLATION .METHOD 
NPTS.UVWPT. ARRAY 
UVWPT. ARRAY 
NPTS.KE. ARRAY 
KE. ARRAY 
NPTS.S. ARRAY 
S. ARRAY 
$END 


«  9*0.0, 

*  0.0, 

=  0.0, 

■  0.0, 

=  0.0, 

*  2,  2, 

*  2,  2, 

=  "LINEAR", 
=  1. 

=  6*0.0, 

=  1, 

=  3*0.0, 

=  1, 

=  0.0, 


$ZONE. GEOMETRY 
ZONE. NUMBER 
CONVERSION. FACTOR 
NUMBER . OF . POINTS . BOTTOM 
INTERPOLATION .METHOD. BOTTOM 

X.  ORIGIN. BOTTOM 

Y.  ORIGIN. TOP 

X.  BOTTOM 

Y.  BOTTOM 

NUMBER . OF . POINTS . TOP 
INTERPOLATION .METHOD . TOP 

X.  ORIGIN. TOP 

Y.  ORIGIN. TOP 
X  .TOP 
Y.TOP 

$END 


1. 

1.0, 

2, 

"LINEAR" , 
0.0, 

0.0, 

50*0.0, 

50*0.0, 

2. 

"LINEAR" , 
0.0, 

0.0, 

50*0.0, 

50*0.0, 


IZONE.MESH 
ZONE. NUMBER 
NUMBER. OF. CELLS. I 
NUMBER. OF. CELLS. J 
ADAPT. ZONE 
X. CENTER 
DELTA. X 

X. STRETCH. MODE 
NUMBER . OF . CELLS . LEFT 
NUMBER . OF . CELLS . LEFTCENTER 
NUMBER . OF . CELLS . RIGHTCENTER 
NUMBER . OF . CELLS . RIGHT 
STRETCH. LENGTH. LEFT 


1. 

21, 

21, 

"NO", 

0.0, 

1.0, 

"STRETCH. FACTORS", 
0. 

0. 

21, 

0, 

1.0, 
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STRETCH . LENGTH . LEFTCENTER 
STRETCH . LENGTH . RIGHTCENTER 
STRETCH . LENGTH . RIGHT 
J. STRETCH. MODE 
NUMBER. OF. J. BLOCKS 
J . BLOCK . NUMBER . OF . CELLS 
J . BLOCK . NUMBER . OF . I . COLUMNS 
J. BLOCK. I. COLUMNS 
J. BLOCK. DELY. BOTTOM 
J. BLOCK. DELY. TOP 
J .BLOCK . STRETCH . FACTORS 
SEND 


1.0, 

1.0, 

1.0, 

"STRETCH. FACTORS", 
1. 

21, 

1. 

1. 

0.0, 

0.0, 

1.0, 


SZONE. OUTPUT 

ZONE. NUMBER  *  1, 

NUMBER. OF. I. COLUMNS  =  0, 

I. COLUMNS  =  10*2, 

NUMBER. OF. TAPS  =  0, 

TAP. I. INDEX  *  20*1, 

TAP. J. INDEX  =  20*1, 

SEND 
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6  DUMP  COMBUSTOR  CALCULATIONS 


In  this  section  the  results  of  several  calculations  using  DUMPSTER  are  described. 
The  first  section  describes  the  results  of  our  study  on  the  use  of  several  flux  functions. 
The  next  section  presents  the  results  of  a  set  of  calculations  of  a  dump  combustor 
with  comparisons  to  experimental  data  and  another  numerical  flow  analysis  program. 
In  the  last  section  the  results  of  the  application  of  DUMPSTER  to  a  dump  combustor 
that  is  closely  coupled  with  a  CD  nozzle  are  presented. 


6.1  Comparison  of  Flux  Functions 


The  LGS  method  requires  a  flux  function  that  is  sufficiently  dissipative  so  that  the 
Gauss-Seidel  iteration  process  converges.  The  LGS  method  was  used  with  the  nat¬ 
urally  dissipative  Steger- Warming  flux  function  and  with  Roe’s  flux  function  with 
eigenvalue  smoothing.  The  added  eigenvalue  smoothing  not  only  helps  the  Gauss- 
Seidel  iteration  converge,  but  also  helps  prevent  aphysical  phenomena  like  shock 
expansions.  The  use  of  added  dissipation  and  the  naturally  dissipative  flux  func¬ 
tion  by  Steger- Warming  raised  concerns  about  the  accuracy  of  the  code  especially  for 
resolving  viscous  shear  layers. 

Additional  numerical  dissipation  tends  to  diffuse  shear  layers.  In  this  study  several 
flux  functions  were  used  to  calculated  a  laminar  boundary-layer  flow  on  a  flat  plate 
for  a  freestream  Mach  number  of  2.0,  laminar  Prandtl  number  of  0.7,  and  a  viscosity 
relation  as  follows. 


This  case  was  selected  because  an  approximate  analytical  solution  is  available  for 
comparison[41].  The  flow  field  schematic  is  shown  in  Figure  10.  A  20x20  mesh 
was  used  with  the  inflow  boundary  of  the  domain  at  the  leading  edge  of  the  plate. 
The  calculated  velocity  and  temperature  profiles  are  compared  to  the  approximate 
analytical  solution  in  Figures  11  through  14.  The  second-order  Roe’s  scheme  with 
k  =  0.0  is  seen  to  give  the  best  agreement  with  the  analytical  solution,  and  the  first- 
order  Steger  Warming  flux  splitting  scheme  the  worst.  The  Steger  Warming  scheme 
is  particularly  poor  in  resolving  the  temperature  profiles  in  all  cases.  The  first-order 
Roe’s  scheme  is  substantially  better  than  the  first  and  second-order  Steger  Warming 
schemes. 

Since  the  LU-SGS  solution  method  does  not  require  a  flux  function  with  as  much 
dissipation  as  the  LGS  method,  the  Steger- Warming  flux  function  is  used  only  for 
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Figure  10:  Schematic  Diagram  of  Test  Case  of  Flat  Plate  Boundary  Layer  with 
Freestream  Mach  Number  of  2.0. 


Figure  11:  Comparison  of  Analytical  and  Calculated  Velocity  Profiles  Using  Steger- 
Warming  Fluxes  for  Flat  Plate  Boundary  Layer. 
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Figure  12:  Comparison  of  Analytical  and  Calculated  Velocity  Profiles  Using  Roe’s 
Fluxes  for  Flat  Plate  Boundary  Layer. 


Figure  13:  Comparison  of  Analytical  and  Calculated  Temperature  Profiles  Using 
Steger- Warming’s  Fluxes  for  Flat  Plait*  Boundary  Layer. 


<)2 


Figure  14:  Comparison  of  Analytical  and  Calculated  Temperature  Profiles  Using 
Roe’s  Fluxes  for  Flat  Plate  Boundary  Layer. 


“roughing  out”  poor  initial  conditions.  For  low-speed  flows  the  LU-SGS  method 
can  be  used  with  either  the  Roe’s  or  Harten-Yee’s  flux  functions  without  significant 
additive  dissipation. 


6.2  Dump  Combustor  Calculations 


In  this  section  calculated  results  from  DUMPSTER  are  presented  for  the  axisymmet- 
ric  non-reacting  dump  combustor  shown  in  Figure  15.  This  configuration  is  a  simple 
50.8mm  radius  round  pipe  dumping  into  a  larger  concentric  pipe  with  a  radius  of 
76.2mm.  The  step  height,  H,  is  the  difference  in  radii  of  the  two  pipes  25.4mm 
(1.0  inch).  Three  different  swirl  numbers  are  considered:  S  =  0.0  (no  swirl),  S  =  0.3, 
and  S  =  0.5.  Calculations  used  the  k-t  turbulence  model  with  the  standard  constants 
and  with  modified  constants.  The  computed  results  are  compared  to  experiment  and 
to  calculations  using  a  different  flow  analysis  computer  code[42]. 


93 


LAST  MEASUREMENT 
STATION 


Figure  15:  Axisymmetric  Dump  Combustor  Geometry  With  Swirler  (bottom)  and 
Without  Swirler  (top). 
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Figure  16:  Two-Zone  Computational  Mesh  for  Dump  Combustor  with  5  =  0. 
6.2.1  No  Swirl:  5  =  0.0 

The  computational  mesh  used  in  this  calculation  is  shown  in  Figure  16.  In  this 
calculation,  a  two- zone  mesh  was  used  consisting  of  42x22  mesh  cells  in  the  lower  zone 
(No.  1),  and  32x12  in  the  upper  zone  (No.  2)  mesh.  The  domain  extends  sixteen  step 
heights  downstream  of  the  step  (the  dump).  At  the  inlet,  uniform  profiles  of  velocity 
direction,  total  pressure,  and  total  temperature  are  imposed  as  boundary  conditions. 
The  upper  and  lower  wall  boundary  conditions  are  free-slip,  non-permeable  walls.  At 
the  subsonic  outflow  boundary  the  pressure  was  fixed  to  F*ex^. 

PT  =  101500.0  Pa 
Tt  =  300.6  °K 

P(>xit  =  101325.0  /^(on  upper  wall) 
u/|V'|  =  1.0 
v/\V\  =  0.0 
w/\V\  -  0.0 

k  —  0.5  m2/s2 
r  =  23.0  m2/s3 


The  Reynolds  number  based  upon  the  combustor  inlet  diameter  (101.6  mm)  was 
approximately  125,000.  Calculated  velocity  vectors,  streamlines,  and  contours  of  tur¬ 
bulent  kinetic  energy  are  shown  in  Figures  17,  18,  and  19.  The  corner  recirculation 
zone  is  about  six  step  heights  long  as  compared  to  eight  step  heights  found  exper- 
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Figure  17:  Calculated  Axial-Radial  Velocity  Vectors  (U-V)  for  S  —  0. 
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Figure  18:  Calculated  U-V  Streamlines  for  5  =  0. 


mentally.  This  ,s  apparently  due  to  the  use  of  the  standard  constants  in  the 
odel  used  m  the  calculation.  Some  numerical  experiments  are  described  in  the  next 
sect, on  that  show  the  sensivity  of  the  corner  recirculation  zone  to  the  k  -,  constants. 


igure  20  shows  profiles  of  measured  and  computed  mean  axial  velocities  at  eight 
axial  locations  UK  is  a  reference  velocity  approximately  560  mm  upstream  of  the 
dump  plane.  U„  has  a  value  of  approximately  19.2  m/s.  Due  to  the  free-slip  boundary 
conditions  used  on  the  combustor  wall,  some  disagreement  is  shown  at  regions  near  the 
dump  plane  where  a  secondary  recirculating  region  is  not  predicted  in  the  computed 
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Figure  19:  Calculated  Turbulent  Kinetic  Energy  Contours  for  5  =  0. 

Figure  21  shows  profiles  of  measured  and  calculated  turbulent  kinetic  energy  ( k )  at 
the  same  axial  locations.  (//*  is  used  to  nondimensionalize  the  kinetic  energy.  From 
these  figures  the  standard  /c-e-model  appears  to  overpredict  the  turbulent  kinetic 
energy  which  results  in  higher  eddy  viscosity  through  out  the  region. 

This  viscosity  overprediction  is  evidenced  by  a  shorter  calculated  corner  recirculation 
region  (Figure  18).  However,  the  overall  physics  of  the  no-swirl  dump  combustor  is 
captured  with  fair  agreement  with  experimental  data  using  the  k-e  model  with  the 
standard  values  for  the  k-e  constants. 

The  deviation  of  the  calculated  u  and  k  profiles  at  X/H  =  0.38  indicated  that  the 
boundary  condition  inflow  velocity  profile  was  not  uniform  as  was  used  in  the  calcula¬ 
tions  above.  To  investigate  the  influence  of  the  upstream  boundary  conditions  upon 
the  flow  this  case  was  rerun  with  the  following  inflow  total  pressure  distribution. 


Radial  Position 
(mm) 

Total  Pressure 

(*»)  B 

00.0 

1.01470x10s 

10.2 

1.01470x10s 

20.0 

1.01465x10s 

25.4 

1.01462x10s 

35.6 

1.01450x10s 

45.7 

1.01405x10s 

48.3 

1.01380x10s 

50.8 

1.01330x10s 

The  above  total  pressure  profile  was  derived  by  trial  and  error.  It  was  adjusted  until 
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Figure  20:  Axial  Velocity  Profile  for  S  =  0.0. 
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the  calculated  axial  velocity  profile  best  fit  the  experimentally  measured  profile  at 
X/H  =  0.38.  The  comparisons  of  the  calculated  profiles  to  the  measured  profiles 
are  shown  in  Figures  22  and  23.  Although  the  axial  velocity  profiles  agree  better  at 
x/H  =  0.38  and  the  turbulenct  kinetic  energy  was  generally  better  predicted,  the 
calculated  axial  velocity  downstream  of  x/H  =  2.0  agreed  less  well.  The  reason  for 
the  deviation  between  the  calculated  and  measured  profiles  is  due  to  a  variety  of 
things  including  at  least:  the  inflow  boundary  conditions,  wall  boundary  conditions, 
turbulence  model,  and  numerical  effects.  To  reduce  the  number  of  variables  in  the 
study  all  following  calculations  were  performed  with  a  uniform  inflow  total  pressure 
profile.  The  fact  that  other  factors  have  significant  influence  upon  the  calculated  flow 
field  must  be  considered  when  comparing  the  calculations  with  experimental  dat.;\ 


6.2.2  Moderate  Swirl:  S  =  0.3 


The  dump  combustor  flow  was  calculated  for  a  swirl  number  of  0.3.  In  this  calculation 
three  zones  were  used  to  construct  the  computational  mesh  shown  in  Figure  24.  The 
effects  of  the  hub  ring  were  modeled  by  using  the  first  zone  (bottom)  (which  consists 
of  42x6  mesh  cells)  in  the  9  mm-high  region  between  the  axis  of  symmetry  and  the 
outer  radius  of  the  hub.  The  mesh  in  zone  2  consists  of  42x19  mesh  cells  and  the  top 
zone  (zone  3)  consists  of  32x12  mesh  cells.  All  zones  extended  approximately  sixteen 
step  heights  downstream  of  the  dump  plane. 

The  inflow  and  outflow  boundary  conditions  were  the  same  as  for  the  S  =  0.0  case 
with  uniform  inflow  boundary  conditions  except  for  the  inflow  swirl  angle  which  was 
set  to  a  swirl  number  of  0.3  between  the  hub  and  the  outer  radius  of  the  inlet.  Free- si  ip 
walls  were  used  on  all  solid  boundaries.  The  direction  cosines  for  the  inflow  velocity 
are  as  follows: 


u/|V?| 

=  0.9164 

v/\V\ 

=  0.0 

w/\V\ 

=  0.4 

Results  for  S  =  0.3  case  are  shown  in  Figure  25  through  28.  Figures  29  and  30 
show  profiles  of  measured  and  calculated  mean  axial  velocities  (u)  and  swirl  veloci¬ 
ties  (to)  at  eight  axial  locations  for  S  =  0.3.  The  agreement  between  the  computation 
and  experiment  is  reasonably  good.  Significant  discrepancies  between  calculated  and 
measured  velocity  profiles  appear  downstream  of  x/H  =  6.0.  These  disrepancies  are 
explained  by  examining  the  turbulent  kinetic  energy  profiles  at  the  same  locations 
(Fig.  31).  The  measured  high  level  of  turbulent  kinetic  energy  near  the  axis  of  sym¬ 
metry  is  not  predicted  by  the  k-t  model  for  downstream  stations.  This  behavior  has 
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Figure  22:  Axial  Velocity  Profiles  for  S  =  0.0. 
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Figure  24:  Computational  Mesh  for  S  =  0.0. 
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Figure  25:  Calculated  Axial-Radial  Velocity  Vectors  (U-V)  for  S  =  0.3. 
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Figure  28:  Calculated  Turbulent  Kinetic  Energy  Contours  for  S  =  0.3. 


also  been  observed  in  the  other  studies  that  were  mentioned  previously.  The  primary 
corner  recirculating  region  length  is  calculated  to  be  about  4.5  step  heights  which  is 
about  the  experimental  measured  value  of  4.3. 

As  in  the  experiment,  no  vortex  breakdown  is  observed  in  the  computed  results.  Even 
though  the  k-t  model  did  not  produce  the  high  turbulent  kinetic  energy  due  to  the 
swirl  near  the  centerline,  the  code  again  simulated  the  general  features  of  the  flow 
reasonably  well. 

The  above  calculations  can  be  compared  to  the  calculations  of  Vanka,  et  al.[42]  shown 
in  Figures  32,  33,  and  34.  Comparison  of  the  two  calculations  is  complicated  because 
of  the  numerous  differences  in  the  computer  codes  and  in  their  application.  One  major 
difference  is  that  in  the  calculations  by  Vanka  et  al.  the  upstream  boundary  conditions 
were  set  at  x/H  —  0.38  using  the  measured  conditions.  Even  with  this  “exact”  inflow 
boundary  condition  the  calculations  are  quite  similar.  Both  calculations  significantly 
underpredict  the  turbulence  kinetic  energy  near  the  centerline  downstream  oi  x/H  = 
1.0.  Computed  turbulent  kinetic  energy  profile  shapes  are  substantially  different 
from  the  measured  profiles  for  both  computer  codes.  In  both  calculations  the  axial 
velocity  profiles  had  a  bump  at  r/H  =  1.4  at  the  downstream  station  ( x/H  —  12). 
In  both  calculations  the  swirl  velocity  in  the  corner  recirculation  region  was  under 
predicted,  probably  as  a  result  of  not  resolving  the  secondary  corner  recirculation. 
Both  calculations  did  very  well  in  simulating  the  swirl  velocity  profiles  upto  x/H  = 
8.0.  Downstream  of  this  station  the  experimental  profiles  of  swirl  were  much  flatter 
than  the  either  of  the  numerical  simulations.  The  strong  similarity  between  the  two 
flow  codes  gives  some  credibility  to  the  present  calculations. 
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Figure  29:  Axial  Velocity  Profiles  for  S  =  0.3. 
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Figure  30:  Swirl  Velocity  Profiles  for  S  =  0.3. 


107 


k/U 


Figure  31:  Turbulent  Kinetic  Energy  Profiles  S  =  0.3. 
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Figure  32:  Axial  Velocity  Profile  Comparisons  S  =  0.3  Calculated  by  Vanka,  et  al.. 
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Figure  34:  Turbulent  Kinetic  Energy  Profiles  S  =  0.3  Calculated  by  Vanka,  et  al.. 
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Figure  35:  Calculated  Axial-Radial  Velocity  Vectors  (U-V)  for  S  —  0.5. 
6.2.3  High  Swirl:  S  =  0.5 


This  case  is  the  the  same  as  the  previous  case  (5  =  0.3)  except  for  that  the  swirl 
number  is  now  set  to  0.5.  The  computational  domain  and  boundary  conditions  are 
the  same  except  for  the  inflow  velocity  angles  listed  below. 


u/\V\ 

=  0.8084 

v/\V\ 

=  0.0 

w/\V\ 

=  0.589 

Results  for  this  calculation  are  shown  in  Figure  35  through  38.  The  length  of  the 
corner  recirculation  zone  is  about  3.5  step  heights  which  is  approximately  equal  to 
the  experimental  value  of  3.2.  A  CTRZ  is  apparent  in  the  exerimental  data  for  this 
case,  but  did  not  appear  in  the  calculation.  However,  the  velocity  profiles  indicate  that 
a  CTRZ  is  close  to  forming.  A  small  recirculation  bubble  is  shown  just  downstream 
of  the  hub.  Figures  39  and  40  show  profiles  of  measured  and  computed  mean  axial 
velocities  (u)  and  swirl  velocities  ( w )  at  eight  axial  locations.  Like  the  5  =  0.3  rase, 
the  k-t  model  under  predicts  the  turbulent  kinetic  energy  near  the  center  line  of  the 
dump  combustor  (Figure  41). 


6.2.4  Calculations  with  Modified  k-t  Constants  for  S  =  0.0. 

Figures  42  through  46  show  the  calculated  results  for  the  no-swirl  case  with  slightly 
modified  k-t  constants.  The  modified  constants  are:  =  0.07  and  Cc i  =  1.887. 
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Figure  38:  Calculated  Turbulent  Kinetic  Energy  Contours  for  S  =  0.5. 

Lowering  the  values  of  these  constants  result  in  lower  values  of  the  eddy  viscosity. 
The  calculated  streamlines  shown  in  Figure  43  show  the  primary  corner  recirculating 
region  is  predicted  very  well,  about  7.5  step  heights  compared  to  8  step  heights  ob¬ 
served  experimetally.  This  calculation  shows  that  the  k-t  turbulence  model  constants 
have  a  significant  influence  on  the  recirculating  flow.  Small  changes  in  the  turbulence 
model  result  in  significant  changes  in  the  flow  field. 


6.2.5  Calculations  with  Modified  k-t  Constants  for  S  =  0.5. 


The  S  =  0.5  case  computed  with  the  standard  k-t  constants  shows  the  largest  dis¬ 
crepancy  with  the  experimental  data,  since  the  computed  results  show  no  vortex 
breakdown.  This  case  was  recalculated  with  modified  turbulence  model  constants: 

=  0.07  and  Ce2  =  1-90.  The  calculated  results  are  shown  in  Figures  47  through  50. 
The  length  of  the  CRZ  region  agrees  very  well  with  the  experimental  results.  The 
extent,  location,  and  size  of  the  CRTZ  is  also  very  well  simulated.  It  is  important  to 
note  that  the  CTRZ  extends  upsteam  of  the  dump  plane.  Comparison  of  calculated 
velocity  profiles  (axial  and  swirl)  and  turbulent  kinetic  energy  profiles  with  measured 
values  in  Figures  51  through  53  show  that  the  calculations  with  the  modified  turbu¬ 
lence  model  constants  agree  favorably  with  experimental  results. 
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Figure  39:  Axial  Velocity  Profiles  for  S  =  0.5. 
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Figure  40:  Swirl  Velocity  Profiles  for  S  =  0.5. 
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Figure  41:  Turbulent  Kinetic  Energy  Profiles  for  S  =  0.5 
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Figure  44:  Calculated  Turbulent  Kinetic  Energy  Contours  with  =  0.07  and  Ca  = 
1.887  for  S  =  0.0. 
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Figure  45:  Axial  Velocity  Profiles.  Calculations  with  CM  =  0.07  and  C(2  =  1.887  for 

S  =  0.0. 
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Figure  46:  Turbulent  Kinetic  Energy  Profiles.  Calculations  with  =  0.07  and 
C'2  =  1.887  for  S  =  0.0. 
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Figure  47:  Calculated  Axial-Radial  Velocity  Vectors  (U-V)  with  Ca  =  0.07  and  Ct 2  — 
1.90  for  S  =  0.5. 
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Figure  48:  Calculated  Swirl-Radial  Velocity  Vectors  (W-V)  with  Cfl  =  0.07  and 
Ct2  =  1.90  for  S  =  0.5. 
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Figure  51:  Axial  Velocity  Profiles.  Calculations  with  CM  =  0.07  and  Cc 2  =  1-90  for 
S  =  0.5. 
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Figure  52:  Swirl  Velocity  Profiles.  Calculatic 


S  =  0.5. 
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Figure  53:  Turbulent  Kinetic  Energy  Profiles.  Calculations  with  =  0.07  and 
C t7  =  1.90  for  S  =  0.5. 
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6.3  Coupled  Dump  Combustor/CD-Nozzle  Calculations 


To  demonstrate  the  capability  of  DUMPSTER  to  calculate  the  flow  in  a  closely  cou¬ 
pled  dump  combustor  and  CD  nozzle  and  to  study  the  influence  of  the  CD  nozzle 
on  the  combustor  flow  field,  four  calculations  were  performed  on  the  configurations 
shown  in  Figure  54.  The  only  differences  between  the  calculations  were  the  lengths  of 
the  combustor  region.  The  combustor /nozzle  geometry  is  similar  to  the  cases  above 
in  that  the  inlet  and  step  height  are  the  same.  Three-zone  meshes  that  were  used 
are  shown  in  Figure  54.  The  dimensions  of  each  mesh  zone  is  shown  in  the  table 
below.  The  combustor  length  (L)  is  the  distance  in  meters  from  the  dump  plane  to 
the  beginning  of  the  nozzle. 


Case  No: 

1 

2 

3 

4 

Combustor 
Length  (L) 

0.1m 

0.2m 

0.3m 

0.4m 

Zone  3 

50x10 

65x10 

80x10 

85x10 

Zone  2 

60x17 

75x17 

90x17 

95x17 

Zone  1 

60x4 

75x4 

90x4 

95x4 

The  inflow  boundary  conditions  are  listed  below  for  an  inlet  swirl  number  of  0.5. 

PT  -  101500.0  Pa 
Tt  =  300.6  °K 
ul\V\  =  0.8084 
v/\V\  =  0.0 
w  l\V\  =  0.589 
k  =  0.5  m2 /s2 
t  =  23.0  m2/s3 


The  outflow  boundary  condition  is  supersonic  outflow.  The  k-t- turbulence  model 
constants  were  set  to  the  modified  values  of  the  last  case. 

The  results  of  the  calculations  are  presented  in  Figures  55  through  59.  Comparison  of 
these  four  cases  illustrate  the  effects  of  the  nozzle  position  upon  the  combustor  flow 
field.  As  the  combustor  region  is  shortened  the  nozzle  effects  the  size  and  extent  of 
the  CTRZ  and  the  CRZ.  As  the  combustor  length  is  decreased  the  CD  nozzle  has  a 
stronger  influence  upon  the  combustor  flowfield.  The  CRZ  increases  in  length  as  the 
combustor  length  is  shortened,  and  i,he  CTRZ  reduces  in  length  and  finally  vanishes. 
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The  CTRZ  changes  continuously  with  no  sudden  change  in  shape  as  the  nozzle  is 
moved  closer  to  the  dump  plane. 

Figures  60  to  62  show  the  calculated  profiles  of  axial  and  swirl  velocities  and  the 
turbulent  kinetic  energy  at  four  stations  along  the  axis  for  all  four  combustor  lengths 
(L).  The  calculations  were  run  3000  steps  with  the  Steger- Warming  flux  function  to 
“rough  out”  the  initial  conditions,  converging  the  solution  three  orders  of  magnitude. 
Then  each  case  was  restarted  using  the  Harten-Yee  flux  function  and  run  6000  steps 
converging  another  three  orders  of  magnitude.  The  overall  run  time  for  these  cases 
was  approximatly  8  CPU  hours  on  an  IBM  RS/6000  model  320  workstation. 
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Figure  56:  Calculated  Swirl- Radial  Velocity  Vectors  (W-V)  for  Dump  Combustor  and 
Closely  Coupled  CD  Nozzle. 
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Figure  57:  Calculated  l 
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Figure  59:  Calculated  Turbulent  Kinetic  Energy  Contours  for  Dump  Combustor  and 
Closely  Coupled  CD  Nozzle.  S  =  0.5. 
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Figure  61:  Swirl  Velocity  Profiles  for  Dump  Combustor  and  Closely  Coupled  CD 
Nozzle.  5  =  0.5 
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CONCLUSIONS  AND  RECOMMENDATIONS 


1.  The  DUMPSTER  computer  code  was  developed  and  applied  to  nonreacting 
axisymmetric  dump  combustor  flow  field  with  closely  coupled  CD  nozzle. 

2.  The  code  uses  a  new  implicit  solution  (LU-SGS)  procedure  which  allows  the 
use  of  low-dissipation  flux  functions  by  Harten-Yee  and  Roe.  Low  dissipation 
is  advantageous,  because  it  introduces  less  aphysical  diffusion  which  can  reduce 
the  accuracy  of  the  simulation.  The  LGS  solution  procedure  was  also  developed. 
The  LGS  was  dependent  upon  a  dissipative  flux  function  for  convergence,  and 
therefore,  discarded  in  favor  of  the  LU-SGS  method. 

3.  DUMPSTER  has  been  used  to  simulate  the  flow  field  of  several  a  non- reacting 
dump  combustor  flow  field.  The  computed  results  agree  resonably  with  both 
the  measurements  and  with  another  flow  analysis  that  had  been  applied  to  the 
same  problem. 

4.  DUMPSTER  has  been  documented  and  delivered  in  source  code  form  for  use  in 
the  AFAPL.  The  code  uses  an  extended  NAMELIST  input  format  to  simplify 
setting  up  problems.  The  code  tracks  the  transport  of  trace  species  and  has 
options  for  both  the  Baldwin-Lomax  algebraic  turbulence  model  and  the  k-t 
two-equation  model.  The  computational  mesh  is  constructed  of  one  or  more 
zones  (or  blocks)  of  mesh  cells,  making  mesh  generation  and  specification  easy. 

5.  The  calculations  indicate  that  the  standard  k-t  turbulence  model  is  not  ade¬ 
quate  for  accurately  simulating  swirling  dump  combustor  flow  fields.  It  does 
provide  a  resonably  good  simulation  if  the  constants  are  adjusted,  but  this  would 
not  be  satisfactory  for  predicting  dump  combustor  flow  fields. 
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